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The  principal  objective  of  this  investigation  was  to  provide  a  reliable  and  efficient  prototype  software 
for  the  engineering  design  and  analysis  of  multilayered  composite  shells,  capable  of  modeling  linear 
and  nonlinear  behavior  in  three  dimensions.  In  addition,  it  should  be  capable  of  assessing  the  quality 
of  the  solution  and  providing  feedback  on  the  basis  of  which  the  solution  quality  can  be  improved.  Its 
hierarchic  structure  should  allow  the  selection  of  models  of  increasing  complexity  in  an  adaptive  way, 
such  that  the  goals  of  computation  are  satisfied  within  the  required  accuracy  and  with  minimal  effort. 

The  project  addressed  the  investigation  of  a  hierarchic  sequence  of  models  for  multilayered  composite 
shells  and  their  implementation  within  the  framework  of  the  p-version  of  the  finite  element  method. 
In  the  Phase  I  project,  we  investigated  the  use  of  the  hierarchic  models  for  the  analysis  of  bending  of 
laminated  (flat)  plates.  Phase  II  utilized  the  results  of  Phase  I  to  fully  assess  the  problems  associated 
with  bending/membrane  coupling  and  curvature,  and  extend  those  results  to  the  nonlinear  (small- 
strain,  large-deformation)  solution  methods  for  laminated  shells,  where  the  most  significant  techno¬ 
logical  contributions  are  to  be  realized.  The  specific  objectives  addressed  in  this  investigation  were: 

•  Investigation  of  the  problems  associated  with  the  implementation  of  a  hierarchic  sequence  of  models 
for  laminated  shells  within  the  framework  of  the  p-version  of  the  finite  element  method  as  an  exten¬ 
sion  of  the  work  done  during  the  Phase  I  project  for  laminated  plates. 

•  Investigation  of  the  added  complexity  for  incorporating  coupling  between  membrane  and  bending 
terms,  caused  by  non-symmetric  stacking  sequences,  curvature,  etc.,  and  the  numerical  generation 
and  orthogonalization  of  the  transverse  shape  functions  for  shells  as  an  integral  part  of  the  solution 
process. 
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•  Investigation  of  the  use  of  the  error  estimators  proposed  during  the  Phase  I  research,  and  the  alter¬ 
nating  projection  method  as  possible  adaptive  strategies  for  the  automatic  selection  of  models  from 
the  hierarchy  for  a  particular  application. 

•  Based  on  the  experience  acquired  with  the  design  of  the  prototype  software  during  Phase  I  research, 
enhance  the  prototype  software  for  shells  to  assess  the  effectiveness  of  the  hierarchic  models  to  solve 
computationally  intensive  problems. 

The  following  objectives  were  achieved  during  the  Phase  II  research  project: 

Task  1:  Investigation  of  the  implementation  issues  of  the  first  model  for  shells  into  a  prototype  soft¬ 
ware.  The  significance  of  this  activity  was  that  it  allowed  addressing  the  finite  element  implementa¬ 
tion  issues  early  in  the  project  and  make  the  necessary  adjustments  in  a  simpler  setting.  Additionally, 
we  were  in  a  position  to  address  the  problems  of  locking  of  shells  and  evaluate  mesh  designs  which 
were  also  very  important  sub-topic  of  the  investigation. 

Task  2:  The  model  selection,  model  construction,  and  generation  of  transverse  shape  functions  for 
higher  order  models  were  investigated.  This  included  the  basic  work  leading  to  adopt  the  proper  strat¬ 
egy  for  model  selection,  for  the  automatic  generation  of  transverse  shape  functions  and  for  the  adap¬ 
tivity  criteria. 

Task  3:  Incorporation  into  the  prototype  software  of  the  additional  models  and  logic  for  model  selec¬ 
tion.  The  additional  models  can  be  manually  or  automatically  selected  from  the  available  set. 

Task  4:  Debugging  and  testing  of  the  prototype  software,  and  solution  of  benchmark  problems,  to 
demonstrate  the  unique  capabilities  of  the  hierarchic  models. 

2.0  The  first  shell  model 


In  the  first  task  of  the  project  we  addressed  the  following  issues: 

•  The  implementation  of  the  first  hierarchic  model. 

•  The  surface  description  of  the  shell  and  the  corresponding  mapping  techniques. 

•  The  formulation  to  account  for  nonlinearities. 

These  points  are  discussed  in  detail  in  the  following. 

2.1  The  first  hierarchic  model 

As  discussed  in  our  Phase  II  proposal,  the  first  hierarchic  shell  model  is  a  five-field  semi-discretiza¬ 
tion  which  approximates  the  solution  corresponding  to  the  three-dimensional  problem  with  the  mini¬ 
mum  number  of  fields. 
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We  consider  a  curved  shell  element  with  arbitrary  geometry  located  in  the  xyz  Cartesian  coordinate 
space.  Let  r,  s,  w  denote  curvilinear  coordinates  such  that  w  =  0  represents  the  middle  surface  of  the 
shell  (Figure  1).  The  shell  under  consideration  is  composed  of  a  finite  number  of  orthotropic  layers  of 
constant  thickness.  The  thickness  of  each  layer  is  denoted  by  /i,-,  and  the  total  thickness  of  the  shell  is 
h  =  hj ,  where  N  is  the  total  number  of  layers.  We  consider  a  displacement  field  which  approx¬ 

imates  tfie  curvilinear  components  of  the  displacement  vector.  Specifically,  the  first  member  of  the 
hierarchic  sequence  of  models  is: 

ur(r,s,w )  =  u\rQ(r,  s)  +  u\[(r,  s)w 

us(r,s,w)  =  u\sQ(r,  $)  +  K|J(r,  s)w  (i) 

uw(r,s,w )  =  M|J(r,j) 

where  (r,  s)  are  the  curvilinear  coordinates  of  the  middle  surface  of  the  shell,  and  w  is  the  direction  of 
the  normal  to  the  middle  surface.  Figure  1  shows  a  typical  quadrilateral  shell  element  in  three-dimen¬ 
sional  space  (Qfc)  and  the  corresponding  standard  quadrilateral  element  (Qs.;). 


FIGURE  1.  Curvilinear  coordinates  associated  with  the  middle  surface  of  the  shell. 


The  main  advantage  in  approximating  the  displacement  components  in  the  curvilinear  system  is  that 
each  field  can  be  augmented  independently  for  the  higher-order  models.  Also,  given  the  variation  in 
the  material  properties  through  the  thickness  in  the  case  of  laminated  composites,  this  makes  it  possi¬ 
ble  to  utilize  a  unique  set  of  transverse  shape  functions  per  field  as  discussed  in  Section  4A  of  the 
Phase  II  proposal.  The  other  advantage  is  related  to  the  specification  of  the  boundary  conditions 
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(loads  and  constraints).  Working  with  the  natural  coordinates  of  the  shell  surface  simplifies  the  imple¬ 
mentation  of  traditional  constraints,  such  as  simple  support,  clamped,  symmetry,  antisymmetry,  etc., 
and  the  specification  of  distributed  surface  or  edge  tractions,  as  well  as  concentrated  forces. 

One  added  complexity  in  working  with  the  curvilinear  coordinates  is  the  incorporation  of  the  rotation 
matrix  into  the  formulation.  The  rotation  matrix  [/?]  is  needed  to  express  the  relation  between  the  glo¬ 
bal  ( xyz )  and  curvilinear  (rsw)  displacement  components, 

{ (2) 

and  its  terms  are  the  components  of  the  normalized  covariant  basis  vectors  (er ,  es ,  evv  in  Figure  1) 
which  are  computed  from  the  derivatives  of  the  mapping  functions.  Introducing  the  following  simpli¬ 
fied  notation  to  indicate  the  difference  between  global  and  curvilinear  components  of  the  displace¬ 
ments:  j  =  u,  u(rsw)  =  u ,  Eq.  (2)  can  be  rewritten  as: 

{u}  =  [i?]{u}  (3) 

The  shell  middle  surface  can  be  written  in  parametric  form  as: 


x  =  x0(r,j),  y  =  y0(r,s), 


z  =  z0(r,s ) 


where  xq,  y0,  zq  are  smooth  mapping  function.  Therefore  considering  the  position  vector  b  shown  in 
Figure  1,  the  normalized  covariant  basis  vectors  are  defined  as: 


vd-y  ds 


ew  “  er  X  ei 


3b  dxo+  dz0f  3b  dxo+  fyo*  dz0f 

dr  ~  dr  1  +drJ  +drk  ’  ds  ~  ds  1  +  ds  J  +  ds  k 

In  Eq.  (6),  ),  j,  k  denote  the  unit  vector  components  in  the  global  coordinate  system.  The  rotation 
matrix  in  Eq.  (3)  is  obtained  from  the  Cartesian  components  of  the  unit  vectors  in  Eq.  (5)  as: 


/?i  Rq  Rrj  erx  e sx  ewx 
[R]=  R2  R5  R&  =  ery  e,y  ewy 

^3  ^6  ^9  ^ rz  ^ sz  ^ wz 
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Note  that  the  rotation  matrix  is  a  function  of  the  curvilinear  components  ( r,  s)  only.  The  relevance  of 
the  presence  of  the  rotation  matrix  in  the  formulation  becomes  apparent  when  considering  the  bilinear 
form  in  the  expression  of  the  principle  of  virtual  work.  Considering  the  case  of  no  body  forces,  no 
spring  boundary  conditions  and  homogeneous  constraints,  the  principle  of  virtual  work  can  be  stated 
as  follows  (Ref.  [1]): 

“Find  [u]  e  S°  such  that  B(u,v)  =  F(v),  for  all  { v}  e  S°” 

where  S°  is  the  space  of  admissible  functions  satisfying  the  homogeneous  boundary  conditions,  {v} 
are  the  test  functions,  F(v)  is  the  virtual  work  of  the  applied  loads,  and  B(u,v)  is  the  virtual  work  of  the 
internal  stresses: 


B{u,v)  =  j([D]{v})T[Q][D]{u}dV 
V 

F(v)  =  \{v}T{T}dA 


In  Eq.  (8),  {T}  is  the  vector  of  the  applied  tractions  in  the  global  coordinate  system,  [D]  is  a  differen¬ 
tial  operator  in  terms  of  the  global  coordinates  and  [Q\  is  the  material  stiffness  matrix.  [Z)]  and  [ Q\  are 
defined  as: 


T 

[D\  = 


Q\\  2 12  2i3  2i4  2i5  Gifi 

i-ooo  4~4~ 

dx  dz  oy 

2 22  Qn  224  Ql5  226 

o  4-  o  4-  o  4- 

[Q]  = 

233  234  235  236 

dy  dz  dx 

(sym)  fi45  246 

0  0  i-i-l-  0 

255  256 

dz  dy  dx 

266 

Introducing  Eq.  (3)  into  Eq.  (8): 

B(u,v)  =  J([T>][R]{v})r[(2][T>]([R]{M}MV 
F(v)  =  \([R]{v})T{T}dA 
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In  Eq.  (10)  the  differential  operator  [ D ]  acts  on  the  rotation  matrix  which  means  that  the  second  deriv¬ 
atives  of  the  mapping  functions  are  also  needed.  The  curvilinear  components  of  the  displacement  vec¬ 
tor  are  approximated  by  polynomial  functions  of  the  form: 

{u}  =  mia},  m  =  m{b}  do 

in  which  [O]  are  known  functions  of  (r,s,w)  and  {a},  {b}  are  the  amplitudes  of  the  basis  functions. 
The  basis  functions  are  given  as  the  product  of  a  function  of  ( r,  s )  times  a  function  of  w  (see  Eq.  (1))  as 
follows: 

1  0  w  0  0 

[$]  =  0  1  0  w  0  (12) 

0  0  0  0  1 

where  ^(r,s)  are  the  hierarchic  basic  shape  functions  characterized  by  the  polynomial  degree  p  and 
the  mapping  functions,  and  are  given  in  Ref.  [1].  Substituting  Eq.  (11)  into  Eq.  (10),  the  expression  of 
the  principle  of  virtual  work  can  be  witten  in  matrix  form: 

f  \ 

{bf  J([T>][R][0])r[e][D][R][0]JR  {a}=  {bf  {T}dA  (13) 

Vy  )  A 

which  has  to  be  satisfied  for  any  {/;}.  Therefore,  Eq.  (13)  can  be  written  in  compact  form  as: 

[£]{a}  =  { q }  (14) 

where  [AT|  is  the  system  stiffness  matrix  and  {q}  is  the  load  vector.  For  any  element  (e)  in  the  mesh, 
the  stiffness  matrix  and  load  vector  terms  are  given  by: 

Kf  =  J([Z)][R]{0}(.)r[e][Z)][R]{0^\/ 

V  (15) 

qf=  limmfiTjdA 
A 

where  {<!>};  is  a  column  of  the  matrix  in  Eq.  (12).  The  solution  of  the  linear  system  of  equations  repre¬ 
sented  by  Eq.  (14),  are  the  curvilinear  components  of  the  displacement  vector. 

As  mentioned  earlier,  the  coefficients  of  the  rotation  matrix  [R]  are  the  Cartesian  components  of  the 
normalized  covariant  basis  vectors  {eresew}^  which  are  computed  from  the  derivatives  of  the  map- 
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ping  functions  of  the  shell  middle  surface.  In  general,  the  components  of  the  normalized  covariant 
basis  vectors  ( er ,  es)  are  non-orthogonal.  In  our  first  implementation  of  this  model  the  non-orthogonal 
basis  was  used.  After  a  few  model  problems  were  solved  using  the  prototype  software,  it  became 
apparent  that  the  computation  of  the  material  matrix  with  respect  to  a  non-orthogonal  basis  would 
impose  additional  complications,  and  a  change  was  required.  An  orthogonal  basis  was  developed  to 
overcome  this  difficulty.  The  orthogonalization  was  performed  by  recomputing  the  unit  vector  in  one 
direction  on  the  shell  surface  as  the  cross  product  of  the  vector  normal  to  the  surface  ( ew )  and  the 
other  unit  vector  on  the  surface  (er):  es  =  ew  x  er  (see  Figure  1).  With  this  change,  the  computation  of 
the  material  stiffness  matrix  was  simplified  for  both  homogeneous  and  laminated  shells. 

The  material  stiffness  matrix  [ Q ]  needed  to  compute  the  element  stiffness  matrix  in  Eq.  (15)  is  deter¬ 
mined  as  follows:  Let  w  (the  shell  thickness  direction)  be  the  direction  of  the  layup  of  the  laminae. 
The  material  properties  of  each  layer  are  defined  in  the  principal  material  directions  of  the  lamina  (. xy 
z).  Let  the  relation  between  the  global  and  lamina  coordinate  systems  at  a  point  within  the  zth  lamina 
be: 


h  h  h 


m  j  m2  m3  i  y 

n ,  n%  z 
L  1  z  JJ  (0 


The  stress-strain  relation  for  the  ith-lamina  in  the  principal  material  directions  is  given  by: 


C ii  Cn  C13  0  0  0 

C22  C23  0  0  0 

C33  0  0  0 


C44  0  0 

C55  0 


(sym) 


or,  in  short  notation, 


{o}(0  =  [C]m{e}(i) 


while  the  same  relation  in  the  global  coordinate  system  is: 


1(0 1  Jxy 


{<*}«)  =  [2](0{e}<0 
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For  an  orthotropic  material,  the  [C]  matrix  contains  9  independent  stiffness  coefficients.  The  relation 
between  the  strains  in  the  lamina  system  and  the  strains  in  the  global  system  is  given  by  (Ref.  [2]): 


l2 

l22 

¥3 

¥  i 

2 

2 

2 

m{ 

m2 

m3 

m2m3 

m3mj 

2 

2 

2 

n\ 

n2 

n3 

n2n3 

n3nl 

2m2n2 

2m3n^ 

m3n2  +  m2n3 

m3ri\  +  m-yn 

2/jftj 

2  l2n2 

2l3n3 

l3n2  -i-  l2n3 

hn\  +  hn2 

2i1m1 

2  l2m2 

2/3m3 

l3m2  +  l2m3 

l3niy  +  lyin-i 

or,  in  short  notation, 

{£}(«)  =  [#](,■){£}(/)  <19> 

The  strain  energy  density  for  the  ith-lamina,  U(i),  is  an  invariant,  and  therefore  it  is  the  same  regard¬ 
less  of  the  coordinate  system: 

U(i)  ~  =  ^ 

Substituting  Eqs.  (17),  (18)  and  (19)  into  Eq.  (20)  we  get: 

um  =  i{e}f0[fi](0{e}(0  =  t{e}(,;.)[if]J)[C](1.)[//](1,{£}(0 


and  therefore,  from  the  above  equation,  the  material  stiffness  matrix  corresponding  to  the  ith-lamina 
in  the  global  coordinate  system  is  given  by: 

[film  =  [H]S)[C](„[H](„  m 

For  each  lamina  the  material  matrix  is  transformed  from  the  material  coordinate  system  to  the  global 
coordinate  system  using  Eq.  (21).  In  computing  the  element  stiffness  matrix  in  Eq.  (15),  the  numerical 
integration  is  performed  layer-by-layer  in  order  to  include  its  material  properties. 
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2.2  Shell  mapping 

The  quality  of  the  mapping  procedure  has  a  substantial  impact  on  the  quality  of  the  finite  element 
approximation.  In  the  p-version  of  the  finite  element  method  large  elements  are  generally  used,  and 
therefore  accurate  representation  of  surfaces  is  essential  so  that  the  errors  of  discretization  can  be  con¬ 
trolled  by  the  mesh  and  the  polynomial  order  rather  than  by  the  mapping  of  the  elements.  We  have 
investigated  a  mapping  technique  for  shells  in  which  the  surfaces  are  approximated  by  piecewise 
interpolation  polynomials  using  special  collocation  points.  This  method,  developed  at  the  University 
of  Maryland  ,  College  Park  [3],  and  investigated  at  Washington  University  in  St.  Louis  [4],  is  called 
Quasi-Re gional  Mapping,  and  has  been  implemented  in  the  prototype  software  for  shells. 

To  demonstrate  the  quality  of  mapping  obtained  by  this  method,  consider  the  problem  shown  in 
Figure  2,  which  represents  the  canopy  of  a  jet  fighter.  The  canopy  was  created  by  connecting  a  set  of 
elliptical  arcs  by  a  NURBS  (Non  Uniform  Rational  B-Spline)  surface.  Two  meshes,  one  consisting  of 
four  quadrilateral  shell  elements,  the  other  of  six  elements,  were  attached  to  the  surface  as  shown  in 
the  figure.  Visually,  the  mapping  is  able  to  capture  all  essential  features  of  the  underlying  surface. 
Numerical  investigation  of  the  quality  of  the  mapping  approximation  and  its  influence  in  the  data 
extraction  from  the  finite  element  solution  can  be  found  in  Ref.  [4]. 

In  the  case  of  shells,  additional  requirements  on  the  quality  of  the  mapping  procedures  are  imposed  by 
the  need  to  approximate  the  second  derivatives  well.  The  second  derivatives  of  the  mapping  functions 
are  required  in  the  computation  of  the  derivatives  of  the  rotation  matrix  in  Eq.  (7)  to  be  used  in  the 
computation  of  the  stiffness  matrix  of  the  elements.  To  illustrate  this  point,  consider  the  product 
[£)][/?]  in  Eq.  (15).  Given  the  definitions  of  [R]  and  [Z)j  in  Eqs.  (7)  and  (9)  respectively,  the  derivative 
of  Ri  with  respect  to  x,  for  example,  will  be  given  as: 

()R  | })[-  dR  ]  f)y 
dx  dr  dx  ds  dx 


but  the  derivative  of  R\  with  respect  to  the  curvilinear  coordinates  should  be  further  expanded  in 
terms  of  the  mapping  functions  in  Eq.  (4)  as  follows: 


d  Xq  dxQf  dxQ  d  x0  dy0  d  y0  dz0  d  Zp 

dR i  =  derx  =  d7  ~5r  t3r  3r2  +3r  df  +3r  dr\ 

dr  dr  zr3/2 


E  = 
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FIGURE  2.  Example  of  Quasi-Regiona! mapping  for  shells. 


and  similarly  for  other  terms.  The  second  derivative  of  the  mapping  function  appears  explicitly  in  the 
derivative  of  the  rotation  matrix  components.  This  clearly  indicates  that  unless  smooth  mapping  func¬ 
tions  are  used,  errors  will  be  introduced  in  the  formulation  which  are  not  related  to  the  those  intro¬ 
duced  by  the  dimensionally  reduced  model.  These  mapping  procedures  have  been  implemented  and 
tested  in  the  prototype  software. 

2.3  Nonlinear  formulation 

Three  types  of  nonlinearities  were  considered  for  the  five-field  semi-discretization  as  indicated  in  the 
following: 

•  Geometric  nonlinearities :  Small-strain,  large-displacement  problems  in  which  the  action  produced 
by  the  loading  changes  as  the  body  deforms. 

•  Eigenvalue  buckling:  Determination  of  the  bifurcation  buckling  loads  for  shells  subjected  to  arbi¬ 
trary  loading. 
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•  Pre-stress  modal  analysis:  Effects  of  pre-stresses  due  to  arbitrary  loading  on  the  natural  frequency 
of  vibration  of  shells. 

The  eigenvalue  buckling  and  the  pre-stress  modal  analysis  are  linearized  nonlinear  problems  that  can 
be  solved  in  a  two-step  operation.  The  geometric  nonlinear  problem  on  the  other  hand  is  a  fully  non¬ 
linear  problem  in  which  the  number  of  required  iterations  is  problem-dependent. 

2.3.1  Geometric  nonlinearities 

Traditional  finite  element  formulations  of  geometrically  nonlinear  problems  are  based  on  the 
Lagrangian  description  of  the  equilibrium  equations.  Several  algorithms  of  this  formulation  are  avail¬ 
able  in  the  literature.  These  include  the  total  Lagrangian,  the  updated  Lagrangian  and  the  ‘co-rota- 
tional  algorithm’  for  beams  and  shells  (see,  for  example,  Refs.  [6]-[7]-[8]). 

The  discretization  procedures  for  geometrically  nonlinear  problems  has  been  associated  almost  exclu¬ 
sively  with  the  h-version  of  the  finite  element  method.  Although  the  p-version  of  the  finite  element 
method  has  been  addressed  in  several  papers  (Refs.  [9]-[10]),  the  formulation  is  based  on  the 
Lagrangian  description. 

In  our  work,  the  geometrically  nonlinear  problem  is  formulated  using  a  weak  form  of  the  spacial/ 
Eulerian  representation  of  the  equilibrium  equations.  This  approach  has  the  following  advantages 
over  the  Lagrangian  description: 

•  The  equilibrium  equations  are  satisfied  in  the  deformed  configuration. 

•  The  displacement  components  approximated  by  high-order  polynomials  and  the  quasi-regional 
mapping  provide  an  accurate  description  of  the  deformed  configuration  and  makes  control  of  dis¬ 
cretization  errors  possible  in  practice. 

•  The  simulation  of  non-conservative  loads,  such  as  follower  loads,  does  not  require  any  additional 
extensions  in  the  formulation,  nor  does  it  affect  the  symmetry  of  the  resulting  stiffness  matrix. 

•  The  extraction  of  stresses  from  the  finite  element  solution  is  straightforward,  since  in  this  approach 
the  equilibrium  equations  are  expressed  in  terms  of  the  Cauchy  stresses. 

The  spatial/Eulerian  formulation  of  the  geometrically  nonlinear  problem  can  be  summarized  as  fol¬ 
lows  (Ref.  [5]):  The  effect  of  large  displacements  is  accounted  for  by  considering  the  equilibrium 
equations  in  the  deformed  configuration.  The  mapping  functions  are  updated  using  the  displacement 
vector,  and  the  strains  and  stresses  are  computed  with  respect  to  the  deformed  configuration  (using  the 
Almansi  strains  and  the  Cauchy  stresses).  The  resulting  nonlinear  system  is  solved  by  iterative  meth¬ 
ods.  For  the  particular  case  of  shells,  the  formulation  is  described  in  the  following. 

Considering  the  case  of  no  body  forces,  no  spring  boundary  conditions  and  homogeneous  constraints, 
the  principle  of  virtual  work  can  be  stated  as  follows: 

“Find  {u}  e  S°  such  that  B(u,v)  =  F(v),  for  all  { v}  <=  5°” 
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where  S°  is  the  space  of  admissible  functions  satisfying  the  homogeneous  boundary  conditions,  {nj 
are  the  trial  functions,  { v}  are  the  test  functions,  F(v)  is  the  virtual  work  of  the  applied  loads,  and 
B(u,v )  is  the  virtual  work  of  the  internal  stresses: 


B(u,v)  =  J{e(v)}r[e]{e(M)Ml/ 
V 

F(v)  =  \{v}T{T}dA 


which  differs  from  Eq.  (8)  in  that  { e(u^ }  is  the  Almansi  strain  tensor  defined  in  the  global  coordinate 
system  due  to  the  trial  functions  and  {e^}  is  the  linear  strain  tensor  due  to  the  test  functions.  That  is: 

{e(M)}  =  {4u)}-{Ae(t0} 

(23) 

{eW}  =  {41”} 

where  {e0)  are  the  linear  components  of  strain  and  { Ae}  are  the  nonlinear  components  as  shown 
below: 


4 (!/)  =  \(Ui,j+Uj,d= 


A  (M)  1 

Aeij  =  2  Wkj’ 


4<i»  =  j(vw  +  vy,,)=  [DKv} 


k  =  1,2,3 


where  the  repeated  index  indicates  summation,  and  utj  represents  the  derivative  of  with  respect  to 
Xj.  Note  that  only  the  linear  definition  of  strains  is  considered  for  the  trial  function.  Substituting  Eq. 
(23)  into  the  bilinear  form  of  Eq.  (22),  we  have: 

B(u,v)  =  J{e(v)}r[(2]{£o“)  }dV-j{EM}[Q]{Ae(u)}dV  (25 

V  V 

and  considering  the  relationship  between  the  curvilinear  and  global  components  of  the  displacement 
vector  in  Eq.  (3)  and  the  strain  definition  in  Eq.  (24),  we  have: 

B(u,v)  =  J([Z)][/?]{v}f[<2][/)]([^]{M})^-J([£>][^]{v})7'[e]{Ae(“)}JV 

V  ~  V  (26 

F(V)  =  J([fl]{v})r{r}<M 
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Substituting  Eq.  (11)  into  Eq.  (26)  and  rearranging,  the  principle  of  virtual  work  can  be  written  as: 


I  ( 

{b}T  f([£>][/?][0])r[(2][Z)][i?][0]jy  {a}=  {b}T  f([i?][0]f{r}iA  + 

W  )  vA 

J([Z)][/?][<l>])r[!2]{Ae(M)}JV 


which  has  to  be  satisfied  for  any  {b}.  Note  that  the  contribution  of  the  nonlinear  strains  have  been 
moved  to  the  right-hand-side  of  the  equation,  as  a  ‘load  vector’  term.  Therefore,  Eq.  (27)  can  be  writ¬ 
ten  in  compact  form: 


[K]{a]  =  {q}  +  {Aq} 


where: 


Kij  =  |([D][/?]{<D},.)7'[e][£>][/?]{<I>^V 

qt  =  lammfiTydA 

A 

A  qt  =  J([D][^]{0}(.)r[e]{  Az{u)}dV 


Note  that  A qt  depends  on  the  solution,  and  therefore  the  system  can  only  be  solved  by  an  iterative  pro¬ 
cedure.  Denoting  the  Ath-iteration  by  a  superscript  (k),  Eq.  (28)  can  be  witten  as: 

[*k-lhiaW}  =  {Jk-l)}  +  {  A^1}}  (29) 

In  the  computation  of  the  stiffness  matrix,  load  vector  and  contribution  of  nonlinear  strains  in  Eq. 
(29),  the  mapping  updated  by  the  solution  corresponding  to  the  (k- 1)  iteration  is  used.  In  other  words, 
during  the  fcth-iteration,  the  mapping  functions  are  given  by: 

x =  xQ(r,s)  +  u ^  ^(r,5,0) 

yW  =  y0(r,^)  +  M^_1)(r,5,0)  (30) 

z(k)  =  Z0(r,s)  +  uf~  1}(r,^,0) 
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where  the  displacements  in  the  global  coordinates  are  computed  from  the  curvilinear  displacement  at 
the  shell  middle  surface  ( w  =  0)  using  the  rotation  matrix  corresponding  to  the  (k-l)  iteration: 


This  formulation  has  been  implemented  in  the  prototype  software  for  planar  and  three-dimensional 
elasticity  problems.  The  results  of  the  preliminary  investigation  clearly  demonstrated  the  potential  of 
the  proposed  method  for  solving  geometrically  nonlinear  problems  for  shell  structures.  An  example  of 
the  implementation  is  included  in  the  next  section. 

2.3.2  Eigenvalue  buckling 

Another  class  of  problems  in  the  analysis  of  shell  structures  is  eigenvalue  buckling  which  is  a  linear¬ 
ized  form  of  a  geometrically  nonlinear  formulation,  useful  for  estimating  the  limits  of  elastic  stability. 
A  capability  to  perform  eigenvalue  buckling  for  homogeneous  and  laminated  shells  was  developed 
and  implemented.  The  main  points  of  the  formulation  are  outlined  in  the  following. 

The  undeformed  configuration  of  the  shell  is  denoted  by  Q  and  its  boundary  by  dfl  The  infinitesimal 
strain  is  defined  in  terms  of  the  Cartesian  components  of  the  displacements  (m,-,  i=1,  2,  3): 

eij  =  \(ui,j  +  uj,d 

which  is  a  simplification  of  the  Green-Lagrange  strains  defined  by 

1, 

“y  ^y  2  M(x> '  y') 

The  simplification  is  justified  by  the  assumption  that  \uq\  «  1  and  hence  the  product  terms  ua  i  uaj 
are  negligible  in  relation  to  iq The  stress-strain  relationship  is: 

aij  ~  cij  +  cijkizki 

where  o^0  is  a  pre-existing  stress  state,  independent  of  w,-,  and  is  the  tensor  of  the  elastic  moduli 
of  the  material.  An  important  property  of  a,y°  is  that  it  is  in  equilibrium  with  the  corresponding  trac¬ 
tions  Tp  =  Oy°  tij  in  the  sense: 

\\G°ij(vi,j+vj,i')dv  =  J  rfvjdA  for  all  vf  e  E  (Q)  (3i) 

q  an 
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where  dV  and  dA  represent  the  differential  volume  and  differential  area,  respectively,  and  E  (Q)  is 
the  space  of  kinematically  admissible  perturbations. 

When  the  reference  configuration  is  stress-free  (i.e.,  Gy°=  0)  then  the  potential  energy  is  defined  by: 

n(«,)  =  \\cllut,fiudV-  J  T,u,dA 

Q.  dQ. 

The  exact  solution  minimizes  II  on  the  set  of  all  kinematically  admissible  functions  denoted  by 
E(Q).  When  the  reference  configuration  is  not  stress-free  then  the  work  done  by  G^°  due  to  the  non¬ 
linear  strain  terms  may  not  be  negligible.  Therefore  the  potential  energy  expression  is  written  in  the 
following  form: 


no,,)  =  J  T,u,M  <32, 

n  n  do. 

The  second  integral  in  Eq.  (32)  represents  the  work  done  by  the  initial  stresses  due  to  the  nonlinear 
strain  terms.  The  work  done  by  Cy°  due  to  the  linear  strain  terms  is  cancelled  by  the  work  done  by  7}° 
in  the  sense  of  Eq.  (31).  The  discretized  form  of  the  potential  energy  in  Eq.  (32)  is: 

n  =  \{af  [K]{a}+l-{af  [G]{a]  -  {af  [q] 

where  { a }  represents  the  coefficients  of  the  basis  functions;  [^f]  is  the  stiffness  matrix;  [G]  is  called 
the  geometric  stiffness  matrix,  and  { q )  is  the  load  vector.  In  typical  structural  stability  problems,  7}=0, 
Gy°  is  predominantly  compressive,  and  the  objective  is  to  find  the  lowest  scalar  multiplier  of  G,y°, 
denoted  by  X,  and  the  corresponding  nontrivial  displacement  vector  function  such  that 

1  X  o 

n(«,-)  =  2  J  Cijkl Eij £kl dV+  2 1  Gij  ua,  i  “a, ;  ^  ^ 

Q.  Q 

is  minimum.  The  stress  field  Gy°  is  called  the  pre-buckling  stress  state,  and  the  critical  load,  also 
called  the  bifurcation  buckling  load,  is  ^mjnT/ • 

The  discretized  form  of  the  eigenvalue  problem  represented  by  the  minimization  of  Eq.  (33)  is: 

([^]  +  X[G]){a}  =  0  (34) 

where  the  stiffness  matrix  is  computed  in  the  usual  way  (see  Eq.  (15))  and  the  geometric  matrix  is 
determined  from  the  second  integral  in  Eq.  (32)  as  follows: 
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{a}T[G]{a}=\(({D}uf[(50]{D}ux  +  ({D}uy)T[o0\{D}uy  +  ({D}uz)T[(50]{D}uz)dV  05) 

£2 

where  ,  uz  are  the  Cartesian  components  of  the  displacements  which  are  related  to  the  curvilin¬ 
ear  components  through  the  rotation  matrix  [/?],  {D}  is  a  differential  operator  vector  and  [ag]  is  the 
stress  tensor  written  in  matrix  form  as  indicated  below. 


{£>} 


X 

X 

xy 

xz 

X 

T  _ 

yx 

y 

yz 

T 

T 

L  zx 

z_\ 

The  implementation  of  this  algorithm  in  the  prototype  software  requires  two  steps  for  the  computation 
of  the  buckling  load  factor: 

•  A  linear  elastostatic  problem  is  solved  first  for  the  specified  loadings  ( T f)  and  constraints.  The  lin¬ 
ear  solution  is  used  to  compute  the  initial  stress  tensor  c(y°.  The  stress  tensor  o^0  at  each  integra¬ 
tion  point  is  used  to  compute  the  geometric  matrix. 

•  After  the  geometric  matrix  is  available,  the  eigenvalue  problem  represented  by  Eq.  (34)  is  solved  to 
find  the  minimum  buckling  load  factor.  The  critical  load  is  then:  Tcr  =  'km\nT{). 

Several  model  problems  of  homogeneous  and  laminated  shells  have  been  solved  using  the  algorithm 
described  above.  The  results  indicate  that  the  implementation  provides  very  accurate  results  when 
compared  with  data  published  in  the  literature  and  with  the  results  of  fully  three-dimensional  analysis. 
An  unique  feature  of  the  formulation  is  that  it  is  not  tied  to  a  particular  type  of  dimensional  reduction 
but  rather  it  can  be  used  in  conjunction  with  the  hierarchic  family  of  models  and  even  for  fully  tree¬ 
dimensional  models. 

2.3.3  Pre-stress  modal  analysis 

The  formulation  for  elastic  vibration  is  analogous  to  Eq.  (32)  for  eigenvalue  buckling.  Specifically, 
we  seek  to  find  03  and  u(e  E  (£2),  m;.  ^  0  such  that 

n*(M;)  =  \  1  Cijkl£tJEkidv  +\\  VV  iuajdv  -  0)2  j  putu,dV  <36) 

a  q  a 

is  minimum.  The  symbols  03  and  p  in  Eq.  (36)  represent  the  natural  frequency  and  the  mass  density, 
respectively.  The  importance  of  the  stress  field  o,y°  is  clearly  visible  from  Eq.  (36):  If  Gty°  is  predomi- 
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nantly  tensile  then  the  stiffness  is  increased,  whereas  if  Gy°  is  predominantly  compressive  then  the 
stiffness  is  decreased.  If  Gy°  is  a  buckling  stress  then  the  lowest  natural  frequency  is  zero. 

The  discretized  form  of  the  potential  energy  in  Eq.  (36)  is: 

n*  =  ±{af  ([K]  +  [G)){a}-(02{a}T  [M]{a}  (37) 

where,  as  before,  {a}  represents  the  coefficients  of  the  basis  functions;  [/£]  is  the  stiffness  matrix;  [G] 
is  the  geometric  stiffness  matrix,  and  [M]  is  the  mass  matrix.  The  discretized  form  of  the  eigenvalue 
problem  represented  by  the  minimization  of  Eq.  (37)  is: 

([K]  +  [G]-k[M]){a}  =  0  (38) 

where  ?i=co2,  the  stiffness  matrix  [£]  is  computed  from  Eq.  (15),  the  geometric  matrix  [G]  is  deter¬ 
mined  from  Eq.  (35)  and  the  mass  matrix  [M]  is  computed  from: 


{a}T[M]{a}  =  Jp {ur  +u2s  +  uw  )dV  <39> 

The  implementation  of  the  pre-stress  modal  analysis  in  the  prototype  software  requires  two  steps: 

•  A  linear  elastostatic  problem  is  solved  first  for  the  specified  loadings  ( T f)  and  constraints.  The  lin¬ 
ear  solution  is  used  to  compute  the  initial  stress  tensor  G^0.  The  stress  tensor  G^°  at  each  integra¬ 
tion  point  is  used  to  compute  the  geometric  matrix. 

•  The  stiffness  matrix  is  modified  by  the  geometric  matrix  and  the  mass  matrix  is  also  computed  to 
solve  the  eigenvalue  problem  of  Eq.  (38). 

Examples  problems  are  presented  in  the  next  section. 

2.4  Example  problems 

2.4.1  Problem  1 :  Linear  elastostatic  analysis  of  a  cylindrical  shell 

A  cylindrical  shell  clamped  at  one  end  and  loaded  by  a  uniform  distributed  normal  traction  was  ana¬ 
lyzed  using  the  5-field  shell  model  implemented  in  the  prototype  software  and  the  results  were  com¬ 
pared  with  a  3D-solid  finite  element  solution  of  the  same  configuration. 

Figure  3  shows  the  shell  model  consisting  of  one  quadrilateral  shell  element  and  the  contour  plot  of 
the  Uy  displacement  component.  The  shell  has  a  radius  R= 2.0,  a  length  a= 2.0,  a  width  b= 0.5  and  a 
thickness  h= 0.10.  The  normal  load  has  a  magnitude  of  q=\00  and  the  material  is  homogenous  and  iso¬ 
tropic  with  £=10xl06  and  v=0.0.  The  corresponding  3D-solid  model,  consisting  of  one  hexahedral 
solid  element,  and  the  contour  plot  of  the  Uy  displacement  component  are  shown  in  Figure  4. 
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FIGURE  3.  Problem  1  -  Mesh  and  contour  plot  for  the  shell  model. 
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FIGURE  4.  Problem  1  -  Mesh  and  contour  plot  for  the  3D-solid  model. 


The  finite  element  solution  was  obtained  for  a  fixed  finite  element  mesh  and  for  increasing  polyno¬ 
mial  order  ranging  from  p=l  to  6  for  both  models.  The  estimated  error  in  energy  norm  for  each  case  is 
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shown  in  Figure  5  in  tabular  form  and  in  a  log-log  scale  plot,  where  the  horizontal  axis  is  the  log  of 


Cfe  fe*'  •k:  M 

Jbmm n«iiiei,MM  %*•*  im.fhmarnb* 


(a)  Shell  model 


Error  Estimate  (All  Elements),  ID-  SOL,  run  #1  to  #6 


Run  # 

DOF 

Total  Potential  Energy 

Rate  of 
Convergence 

Estimated 
%  Error 

1 

10 

-2 . 389949952653914e-02 

0.00 

99.72 

2 

25 

-4 . 120727750166704e-01 

0.05 

95.05 

3 

40 

-3 . 285389747569158e+00 

1.45 

48.02 

4 

60 

-4 . 259780965826461e+00 

5.63 

4.90 

5 

85 

-4 . 269957708088871e+00 

7.40 

0.37 

6 

115 

-4 . 270015928027199e+00 

7.40 

0.04 

Estimated  Limit 

-4 . 270016597868776e+00 

(b)  3D-solid  model 


Error  Estimate  (All  Elements),  ID=  SOL,  run  #1  to  #6 


Run  # 

DOF 

Total  Potential  Energy 

Rate  of 
Convergence 

Estimated 
%  Error 

1 

12 

-2 . 115019881502798e-02 

0.00 

99.75 

2 

36 

-1 . 379345240360249e-01 

0.01 

98.37 

3 

60 

-2 . 589890171918026e+00 

0.88 

62.76 

4 

99 

-4 . 212006486270239e+00 

3.32 

11.93 

5 

153 

-4 . 270967693800999e+00 

4.03 

2.07 

6 

225 

-4 . 272710128732996e+00 

4.03 

0.44 

Estimated  Limit 

-4 . 272791785275626e+00 

FIGURE  5.  Problem  1  -  Estimated  relative  error  in  energy  norm. 


the  number  of  degrees  of  freedom  (DOF)  and  the  vertical  axis  is  the  log  of  the  percent  estimated  rela¬ 
tive  error  in  energy  norm.  Note  that  the  DOF  for  the  shell  model  is  smaller  than  for  the  3D-solid 
model  for  the  same  run  number  (p-level),  and  the  rate  of  convergence  increases  substantially  for  p- 
level  greater  than  or  equal  to  4.  The  shell  model  is  converging  to  the  same  potential  energy  as  the  ref¬ 
erence  three-dimensional  solution  given  by  the  3D-solid  model. 

The  displacement  of  point  A  (Figure  3)  at  the  free  end  of  the  shell  mid-surface  is  shown  in  Figure  6 
for  the  shell  model  (Figure  6a)  and  for  the  reference  3D-solid  solution  (Figure  6b)  as  a  function  of  the 
run  number.  The  results  for  both  models  are  almost  identical.  The  difference  between  the  shell  model 
and  the  3D-solid  model  is  only  0.04%  for  the  Ux  displacement  component  and  0.05%  for  Uy  compo¬ 
nent.  Figure  7  shows  the  convergence  characteristics  of  the  displacement  component  Uy  as  a  function 
of  the  number  of  degrees  of  freedom  (DOF)  for  the  shell  model. 
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Point  Function,  ID=  SOI,  run  #1  to  6 


(a)  Shell  model 


4,  x=  1.68294e+00,  y=  1.08060e+00,  z=  0.D0000e+00,  Ux=-5 . 46170e-04  ,  Uy=-9 . 18733e-04 


4,  x=  1.68294e+00,  y=  1.08060e+00,  z=  0.00000e+00,  Ux=~l . 49293e-02  ,  Uy=-1 . 93801e-02  | 

4,  x=  1.68294e+00,  y=  1.08060e+00,  z=  0.00000e+00,  Ux=-1 . 06430e-01  ,  Uy=-1 . 4642Se-01  ! 

4,  x=  1.68294e+00,  y=  1.08060e+00,  z=  0.00000e+00,  Ux=-1 . 24992e~01  ,  Uy=-1 . 81S03e-01 

4,  x=  1.68294e+00,  y^=  1.08060e+00,  z=  O.OOOOOe+OO,  Ux=  1 . 24740e-01  ,  Uy=-1 . 81S82e-01  | 

4,  x=  1.68294e+00,  y=  1.08060e+00,  z=  O.OOOOOe+OO,  Ux=-1 . 24741e-01  ,  Uy=-1 . 81B84e-01 { 


Point  Function,  ID=  SOL,  run  #1  to  6 


(b)  3D-Solid  model 


2,  x=  1.68294e*00,  y=  1.08060e+00, 
2,  x=  1.68294e+00,  y=  1.08060e+00, 
2,  X“  1.68294e+00,  y^  1.08060e+00, 
2,  x=  1.68294e+00,  y=  1.08060e+00, 
2,  x=  1.68294e+00,  y=  1.08060e+00, 
2,  x=  1.68294e+00,  y=  1.08060e+00, 


8.7S886e-lS,  Ux=-4 . 900S3e-04 
8 . 7S886e-lB ,  Ux— 4 . 28834e-03 
8.7S886e-lB,  Ux=-9 . 11811e-02 
8.7B886e-lB,  Ux=-1 . 2S888e-01 
8.7B886e-lB,  Ux=-1.24780e-01 


Uy=-7 . 21372e-04 
Uy=-B . 92BlBe-03 
Uy=-1.20099e-01 
Uy^-1 . 80942e-01 
Uy=-1 . 81632e-01 


8.7B886e-lB,  Ux=-1 . 24793 e~01  ,  Uy^-1.81682e-01 


FIGURE  6.  Problem  1  -  End  displacement  components  Ux  and  Uy. 


FIGURE  7.  Problem  1  -  Convergence  plot  of  Uy  for  the  shell  model. 


This  example  demonstrates  some  of  the  key  features  of  the  implementation  of  the  first  shell  model 
into  the  prototype  software:  Quality  of  approximation  of  the  three-dimensional  problem;  global  error 
assessment  capability;  and  local  error  assessment  thorough  convergence  checks. 

2.4.2  Problem  2:  Geometric  nonlinear  analysis  of  a  slab 

Consider  the  case  of  a  rectangular  slab  clamped  along  one  edge  and  loaded  by  a  uniformly  distributed 
normal  traction,  as  shown  in  Figure  8a.  The  load  remains  normal  to  the  plate  surface  during  the  defor- 
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FIGURE  8.  Problem  2  -  Geometric  nonlinear  analysis. 


mation  (‘follower’  load).  Figure  8b  shows  the  deformed  configuration  of  the  plate  (full  scale)  for  the 
linear  solution.  It  is  clear  form  the  results  that  the  deformation  obtained  from  the  linear  solution  is  out¬ 
side  of  the  range  of  validity  of  the  linear  theory  of  elasticity.  Figure  8c  shows  the  deformation  (full 
scale)  obtained  with  the  geometrically  nonlinear  formulation  implemented  in  the  prototype  software. 

2.4.3  Problem  3:  Buckling  and  modal  analysis  of  a  roof  structure 

Consider  the  cylindrical  roof  structure  shown  in  Figure  9.  The  shell  is  of  sandwich  construction,  with 
the  outer  layers  of  high  modulus  graphite/epoxy  composite  material  and  an  isotropic  material  core.  A 
vertical  dead  load  is  applied  to  the  roof,  and  the  objective  of  the  analysis  is  to  determine  the  maximum 
vertical  deflection,  the  largest  normal  stress,  the  first  natural  frequency  of  vibration  and  the  buckling 
load  factor.  Also  of  interest  is  to  determine  the  change  in  natural  frequency  as  a  function  of  the  load 
magnitude. 
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FIGURE  9.  Problem  3  -  Cylindrical  roof  structure.  Notation. 

The  following  material  properties  were  used  for  the  external  layers  and  core: 

Graphite/Epoxy:  EL= 25xl06;  jE^ =lxl06;  GLT=  5xl05;  Gjj  =2xl05;  vLT=  v  77  =0.3,  p=lxl0'4 
Isotropic:  E=3xl06;  v=0.0;  p=lxl0'4 

where  L  indicates  the  direction  parallel  to  the  fibers  and  T  is  the  transverse  direction.  In  this  problem, 
the  L-direction  is  aligned  with  the  global  Z-axis.  Because  of  symmetry,  only  one  fourth  of  the  roof 
was  included  in  the  analysis.  The  problem  was  solved  using  the  first  hierarchic  shell  model  and  also 
using  a  3D-solid  model  in  order  to  have  a  reference  solution.  Figure  10  shows  the  2-element  mesh  for 
the  shell  model  and  the  6-element  mesh  for  the  3D-solid  model.  In  the  solid  model,  each  layer  was 
discretized  using  hexahedral  elements. 

Symmetry  boundary  conditions  (un= 0,  where  un  is  the  displacement  normal  to  the  edge)  were  speci¬ 
fied  along  two  orthogonal  directions,  and  antisymmetry  boundary  conditions  (uw=ut=0,  where  ut  is 
the  displacement  tangent  to  the  edge)  were  used  to  represent  the  effects  of  the  end  support.  The  other 
edge  is  free.  Note  that  a  thin  element  was  defined  along  the  free  edge  of  the  shell  and  3D-solid  models 
to  properly  account  for  boundary  layer  effects.  The  results,  obtained  from  shell  model  1  and  from  the 
3D-solid  model,  are  summarized  in  Table  1.  They  include  the  maximum  vertical  displacement, 
wyA(xA,yA,ZA)>  where  (xA,  yA,  zA)  are  the  coordinates  of  the  shell  mid-surface  at  point  A  (see  Figure  9); 
the  normal  stresses  at  point  A  on  the  lower  surface  of  the  roof,  czA;  the  first  natural  frequency  without 
pre-stress,/!,  and  the  buckling  load  factor,  BLF.  The  last  row  in  the  table  are  the  value  of  the  critical 
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load  on  the  roof,  which  were  computed  as  the  product  of  the  applied  load  (PY=0.625)  and  the  buck¬ 
ling  load  factor  ( BLF ). 


table  i .  Results  for  problem  3 


Function 

Shell  model  1 

3D-solid  model 

Difference  (%) 

WyA 

-2.73 

-2.82 

-3.2 

azA 

10614 

10843 

-2.1 

fl  [Hz] 

4.71 

4.64 

1.5 

BLF 

11.29 

10.39 

8.7 

(py)cr 

7.06 

6.49 

8.7 

The  results  of  shell  model  1  are  very  close  to  those  obtained  using  the  3D-solid  model  specially  for 
the  linear  and  modal  analysis  results.  Figure  11  shows  the  buckling  mode  shape  for  both  models.  The 
effect  of  the  pre-stress  induced  by  the  applied  load  (PY)  on  the  first  natural  frequency  of  the  roof 
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Shell  model  3D-solid  model 


FIGURE  11.  Problem  3  -  Buckling  mode  shape 

structure  was  obtained  for  the  shell  model,  and  the  results  are  included  in  Table  2.  The  pre-stress 


table  2.  Effect  of  pre-stresses  on  the  first  natural  frequency 


Py 

0.0 

1.0 

2.0 

3.0 

4.0 

5.0 

6.0 

6.5 

7.0 

7.06 

/l  [Hz] 

4.71 

5.39 

5.82 

6.01 

5.89 

5.31 

4.07 

3.02 

0.95 

0.00 

increases  the  natural  frequency  until  the  load  is  about  half  the  critical  load,  and  then  decreases  rapidly 
as  the  critical  buckling  load  is  approached. 

3.0  Higher  order  models _ 

Hierarchic  sequence  of  models  satisfy  the  equilibrium  equations  of  three-dimensional  elasticity  to  the 
desired  degree  of  accuracy.  In  the  limit  it  converges  to  the  fully  three-dimensional  solution.  Depend¬ 
ing  on  the  goals  of  computation,  the  analyst  can  select  the  model  that  best  fits  the  goals.  Choosing 
progressively  higher  models,  the  computational  effort  increases,  but  of  course  the  accuracy  in  the 
results  is  improved  also.  If  only  structural  response  is  required,  a  low-order  model  is  generally  suffi¬ 
cient,  specially  for  thin  shells.  Hierarchic  sequences  of  models  make  adaptive  selection  of  the  model 
which  is  best  suited  for  the  purposes  of  a  particular  analysis  possible. 

3.1  Concept  and  terminology 

The  first  rigorous  proof  of  the  relation  between  the  three-dimensional  solution  and  a  plate  model  was 
given  by  Morgenstern,  Ref.  [11],  in  1959.  The  construction  of  hierarchic  models  for  homogeneous 
isotropic  plates  and  shells  was  discussed  by  Szabo  and  Sahrmann  in  1988,  Ref.  [12].  The  optimality 
conditions  for  the  construction  of  the  hierarchic  models  for  homogeneous  plates  was  first  discussed 
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by  Schwab  in  1989,  Ref.  [13].  The  extensions  of  these  concepts  for  laminated  composites  was  first 
addressed  in  Ref.  [14]  for  plates  in  cylindrical  bending,  in  Ref.  [16]  for  general  plates  with  mid-plane 
symmetry,  and  in  Ref.  [17]  for  laminated  shells. 

Hierarchic  modeling  terminology  has  been  used  by  other  authors  (see  for  example  Ref.  [19])  but  in  a 
context  which  is  different  from  ours.  As  we  understand  it,  a  hierarchic  sequence  of  models  is  properly 
formulated  if  satisfies  the  condition  that  the  corresponding  exact  solutions  uEX(HMh\  converge  to  the 
exact  solution  of  the  fully  3D-problem  uEX(3D),  for  a  fixed  laminate  thickness, 


lim 

i  — »  o o 


(3D) 

lEX 


(. HM\i)\ 
~  UEX 


E(Cl) 


=  0 


where  E(Q)  is  the  energy  norm.  In  addition,  they  should  posses  two  highly  desirable  features: 

1 .  The  exact  solution  of  each  model  converges  to  the  same  limit  as  the  exact  solution  of  the  corre¬ 
sponding  3D-problem  with  respect  to  the  laminate  thickness  (h)  approaching  zero: 


lim 

h  — ^  0 


(3D)  (HM\i) 

UEX  ~ UEX 

E(Cl) 

(3D) 

UEX 

E(  Q) 

i  =  1,  2, 


This  requirement  is  important  because,  typically,  uE^3D^  in  the  interior  regions  of  the  domain  behaves 
as  if  h  were  close  to  zero,  assuming  that  the  loading  is  smooth. 

2.  Optimality  of  convergence:  When  the  exact  solution  uEX3D^  is  sufficiently  smooth: 


(3D)  (HM\i) 

UEX  ~ UEX 

E(Q) 

u(3D) 

UEX 

E(Q) 

»  Ch1, 


where  C  is  a  constant,  independent  of  v,  y,  the  rate  of  convergence,  is  a  constant  which  depends  on  i, 
and  yi+1  >  yr 

The  hierarchic  sequence  utilized  for  the  Phase  I  project  for  laminated  plates  and  the  one  developed 
during  Phase  II  for  laminated  shells,  satisfy  these  requirements  explicitly. 

3.2  Transverse  shape  functions  for  hierarchic  shell  models 

Two  high-order  (hierarchic)  shell  models  for  homogeneous  and  laminated  composites  were  developed 
and  implemented  within  the  framework  of  the  finite  element  software  product  Stress  Check.  The 
transverse  shape  functions  developed  for  laminated  plates  (see  Ref.  [17])  were  reevaluated  and  modi¬ 
fied  to  be  used  for  laminated  shells.  Six  new  shape  functions  and  their  derivatives  are  available. 
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As  described  in  Section  2.0,  the  first  hierarchic  shell  model  is  a  five-field  semi-discretization  which 
approximates  the  curvilinear  components  of  the  displacement  vector  as  indicated  in  Eq.  (1).  The 
expansion  of  the  three  displacement  components  to  include  two  additional  models  can  be  written  in 
general  form  as: 

ur(r,  s,  w)  =  u\rQ(r,  ^)^>1(w)  +  u\r{(r,  s)®3(»  +  u\r2(r,s)®6(w)  +  u\r3(r,s)®9(w) 
us(r,s,w)  =  a|p(r,  6-)02(w)  +  M|j(r,  5)O4(w)  +  n|2(r,j')O7(w)  +  M|3(r,5)0)10(w)  <40) 

uw(r,s,w )  =  «|o(r,5)05(w)+  w|^(r,  5)0>8(w)+  wf  0vy)<f>n(w) 

where  the  <E>j(w)  are  the  transverse  shape  functions  (director  functions)  which  depend  on  the  variable 
normal  to  the  shell  middle  surface  and  on  the  type  of  shell  (homogenous  or  laminated).  Note  that  there 
is  a  change  in  notation  with  respect  to  that  shown  in  Eq.  (1).  The  unknown  functions  of  (r,  s)  are  iden¬ 
tified  with  a  subscript  that  indicated  a  sequential  order,  and  by  a  superscript  which  refers  to  the  corre¬ 
sponding  displacement  component,  while  the  known  functions  of  w  are  numbered  to  reflect  the 
sequence  of  model  construction. 

The  transverse  shape  functions  are  obtained  according  the  following  procedure: 

Let  ur(r,  s,  w),  us(r,  s,  w),  uw(r,  s,  w)  denote  the  displacement  field  for  the  shell  composed  of  an 
arbitrary  number  of  orthotropic  layers  bonded  together  subjected  to  normal  surface  loading  7}(r,  s ) 
and  satisfying  the  equilibrium  equations  of  3D-elasticity.  Let  the  stresses  be  related  to  the  strains  by 
the  generalized  Hooke’s  law,  and  the  strains  related  to  the  displacements  by  the  small  strain  theory. 
The  problem  is  to  find  the  displacement  field  that  minimizes  the  potential  energy  functional  Il(n)  over 
the  subspace  En(Q),  defined  as: 

n\  n2 

En(Q)  =  <  u\ur(r,s,w)=^  u\j(r,s)<Pj(w),  us(r,s,w)=YJ  u\ j(r,s)\\fj(w), 

j  =  0  7  =  0 

"3 

uw(r,s,w)= u\J(r,s)pj(w ) 

7  =  0 

The  functions  cpy(w),  \|/.(w),  p;(w)  are  derived  on  the  basis  of  the  degree  to  which  the  equilibrium 
equations  of  3D-elasticity  are  satisfied.  The  procedure  to  obtain  these  functions  is  similar  to  the  one 
outlined  in  Refs.  [14]  -  [18]  for  laminated  plates  in  bending.  The  main  step  in  the  derivation  for  lami¬ 
nated  shells  are  as  follows: 
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•  Perform  a  partial  Fourier  transform,  with  parameters  (3  over  the  domain  Q  of  the  three-dimensional 
problem  described  above.  This  transformation  eliminates  two  field  variables  (r,  s )  in  the  displace¬ 
ment  components,  so  that  derivatives  with  respect  to  r  and  s  become  multiplications  by  ;|3  in  the 
Fourier  transformed  variables. 


cp((3,w)  =  | ur(r,s,w)e~l^('r  +  s^dA 

A 

\|/((3,w)  =  jus(r,s,w)e~l^r+s^dA 

A 

p(P,w)  =  J«w(r>s,w)c-iP(r+,>dA 

A 

•  Write  down  the  strain-displacement  relations  in  the  transformed  variables,  substitute  them  into  the 
stress-strain  relations  and  express  the  equilibrium  equations  in  their  Fourier  form.  A  system  of 
ordinary  differential  equations  is  obtained  in  the  variable  w. 

•  Expand  the  functions  cp((3,w),  \|/([3,vv’),  p((3,w)  in  powers  of  P  around  (3=0. 

<p(J3,w)  =  Cp0(w)  +  Pcp^w)  +  P2(p20)  +  P3(p3(w)  +  ••• 

Y(P,w)  =  v0(w)  +  ftyi(w)  +  P2V2<»  +  P3¥3(w)  +  ••• 

p(p,w)  =  p0(w)  +  pcp1(w)  +  p2p2(w)  +  p3p3(w)+ ... 


•  Replace  the  expanded  functions  into  the  Fourier  form  of  the  equilibrium  equations,  which  must  be 
satisfied  for  any  power  of  (3.  The  transverse  shape  functions  are  obtained  by  solving  these  equa¬ 
tions.  For  example,  the  equilibrium  equations  corresponding  to  (3°  are: 

«245Vo  +  255<Po  =  0 
(Q44Vo  +  Q45%Y  =  0 
(Q33P0)  =  0 

•  After  the  functions  are  obtained  by  integration  of  the  equilibrium  equations,  they  are  normalized  in 
such  a  way  that  all  functions  are  zero  at  w=0  and  also  a  thickness  factor  is  included.  For  consis¬ 
tency  of  notation,  the  normalized  transverse  functions  are  defined  as  0,(w),  i'=l,  2,  3, ... 
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The  transverse  shape  functions  for  laminated  shells  depend  on  the  stacking  sequence  layup  and  on  Q y, 
which  are  the  coefficients  of  the  3D-lamina  material  matrix  in  the  shell  principal  directions. 

For  shell  model  1,  the  transverse  shape  functions  are  the  same  for  laminated  or  homogeneous  mate¬ 
rial.  Referring  to  Eq.  (1),  they  are  given  by: 

^(w)  =  ®2(w)  =  <f>5(w)  =  1,  <J>3(w)  =  <E>4(w)  =  w(h/ 2) 

where  h  is  the  shell  thickness,  and  -1  <  w  <  1.  For  the  higher  order  models,  the  transverse  shape  func¬ 
tions  are  defined  as  follows: 

For  homogeneous  shells 

<E>6(w)  =  <l>7(w)  =  w2(h/ 2)2,  <f>8(w)  =  w(h/2) 

<J>9(w)  =  0>10O)  =  w3{h/ 2)3,  On(w)  =  w2(h/ 2)2 

For  laminated  shells 

$6(w)  =  ^{(p2(w)-(p2(°)}  09(w)  =  cp3(w) 

H  f  h\^ 

^7(w)  =  2^20)-  ¥2(°)}  <I>io(w)  =  UJ  V3(w) 

®8(w)  =  - Pi(0)}  ®n(w)  =  (jf)  p2(w) 

W 

f  \  f  2 55 _  ^45  j 

\|/2(w)  =  I - -dw 

_j  ^44^55  “  @45 

Pi(w)  =  f  - ^-dw 
-1 
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It  is  clear  from  the  above  expressions  that  if  <2,;/  are  constant  within  each  lamina  but  different  from 
lamina  to  lamina,  the  transverse  variation  of  cp2(w),  \|/2(w),  p^w)  is  piecewise  linear,  and  the 
slope  ate  each  interface  depends  on  the  material  properties  of  each  layer. 

The  other  three  functions  are: 


(p3(vtO  =  J 


f  ,  .  e44Ya(M')  245 ybM 

Pl(w)  + 


-i 


644655  “645  644655  -Q45  J 


dw-j 


'  ,  ,  644Ya(w)  645Yi(w) 


v3(w)  =  J 


,  ,  QislbW  645  Ya(w) 

Pl(w)  + 

V 


-1 

0 


Q44Q55-Q45  044^55“  Q45  ) 


dw-  j 


Q44Q55-Q45  Q44Q55-Q45  J 


f  ,  ,  QsstbW  245Yfl(^) 

Pi(w)  + 


\dw 


-1 


044^55  ~  £45  Q44Q55  “  Q45  J 


\dw 


f  2w  (  (2l3  +  236)92W 

(Q23  +  | 

r  (2w  (  (Q\3  + 

(Q23  +  236)¥2(w)>j 

V233  G33 

G33  J  J 

l<233  Q33 

233  J 

dw 


-1 


-l 


2i3  +  236,  r2i3  +  236 


6 


dw, 


33 


(*623 +  636  f  623  + 636  j 

=  j  — qZ — rfw,_  J  — - - dw 


-1 


33 


-1 


633 


In  this  case,  the  transverse  variation  is  piecewise  quadratic. 

Several  procedures  have  been  evaluated  for  selecting  the  optimal  distribution  of  fields  for  a  given 
problem.  Starting  with  the  minimal  number  of  fields  (5-field  model),  the  question  is  how  to  construct 
the  next  model?  Our  research  indicated  that  from  an  implementation  and  performance  point  of  view, 
the  next  model  should  be  constructed  by  adding  one  more  field  to  each  displacement  component  (un 
us,  uw),  and  that  the  transverse  function  for  the  r  and  s  components  should  have  the  same  power  of  w. 
In  other  words,  model  2  should  be  an  8-field  semi-discretization  with  the  selection  of  the  three  new 
fields  from  the  available  set  in  Eq.  (40). 

3.3  Automatic  selection  of  models 

The  optimal  selection  of  a  particular  model  form  the  hierarchic  family  of  models  is  problem-depen- 
dent.  Starting  with  a  5-field  model,  the  next  model  is  constructed  by  adding  one  more  field  to  each 
displacement  component,  in  such  a  way  that  model  2  has  eight  fields  and  model  3  has  eleven.  Several 
methods  were  evaluated  to  select  the  best  shape  functions  from  the  available  set  to  construct  the  next 
higher  order  model:  The  alternating  projection  method;  the  residuals  method  and  the  potential  energy 
methods.  These  methods  are  described  in  detail  in  the  following. 
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3.3.1  The  alternating  projection  method 

The  alternating  projection  consists  of  three  parts:  First,  the  solution  corresponding  to  the  lowest  mem¬ 
ber  of  the  hierarchy  is  computed  (the  5-field  semi-discretization).  Second,  the  part  of  the  solution 
which  depends  on  the  transverse  variable  (w  in  Eq.  (1))  is  expanded  in  terms  of  the  known  transverse 
shape  functions  over  each  element.  A  one-dimensional  problem  is  solved  to  find  which  of  these  terms 
in  the  expansion  may  contribute  to  a  significant  change  in  the  potential  energy  of  the  solution.  Third, 
the  model  order  within  each  element  of  the  mesh  is  adjusted  based  on  the  results  of  Part  2,  and  the 
problem  is  solved  again  for  the  new  model  distribution.  These  three  steps  are  now  described  in  detail: 

Step  1:  Start  with  the  5-field  model  for  all  the  elements  in  the  mesh.  Obtain  the  finite  element  solution 
corresponding  to  this  model  as  described  in  Section  2.0.  Note  that  in  this  step  we  solve  a  two-dimen¬ 
sional  problem,  because  the  transverse  variation  of  displacements  is  know  a-priori.  The  results  of  this 
step  are  the  coefficients  in  Eq.  (14).  Therefore,  for  every  element  in  the  mesh  we  know: 

s),  u\sQ(r,s),  «|J(r,  s),  w|*(r,  s),  s) 

Step  2:  Expand  the  displacement  field  given  in  Eq.  (1)  in  terms  of  the  transverse  shape  functions 
shown  in  Eq.  (40)  over  each  element  of  the  mesh: 

ur(r,  s,  w)  =  M|^(r,  ^)  +  w|'1(r,  ^)(d»3(w)  +  af)06(w)-l-a^)09(w)-l- ...) 

us(r,s,w )  =  m|q(t,  5)  +  «|j(r,  ^)(O4(w)  +  pf)O7(w)  +  |32^)O10(w)+ ...)  (4i) 

uw(r,  s,  w)  =  u\™(r,  s)(®5(w) +  y\k)®8(w)  +  +  •••) 

The  transverse  shape  functions  d>z‘(w)  are  those  obtained  in  Section  3.2.  Therefore,  the  only  unknowns 
in  Eq.  (41)  are  the  coefficients  oq,  P,-  and  y(-  that  multiply  these  functions.  The  superscript  k  indicates 
the  &th  finite  element. 

This  new  displacement  field  is  then  used  in  Eq.  (10)  to  obtain  a  new  system  of  equations  to  solve  for 
the  unknown  coefficients  in  Eq.  (41).  Note  that  the  new  system  corresponds  to  a  one-dimensional 
problem  in  the  transverse  direction  since  the  only  unknowns  in  the  expanded  displacement  field  are 
the  coefficient  that  multiply  the  transverse  functions.  The  system  size  is  3 xNxM,  where  M  is  the  num¬ 
ber  of  elements  in  the  mesh,  and  N  is  the  total  number  of  unknown  coefficients  in  Eq.  (41).  The  rela¬ 
tive  magnitudes  of  the  coefficients  will  indicate  which  additional  terms  will  contribute  most 
significantly  to  the  potential  energy  of  the  solution. 

Step  3:  Expand  the  model  order  over  each  element  or  groups  of  elements  according  to  the  relative 
values  of  the  coefficients  cq,  P(-  and  that  multiply  the  transverse  functions.  Once  the  number  of 
fields  are  known,  the  solution  of  the  corresponding  2D  problem  is  obtained  for  the  selected  higher- 
order  model. 
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The  name  alternating  projection  comes  from  the  fact  that  the  three-dimensional  problem  is  solved  by 
first  minimizing  the  potential  energy  of  a  two-dimensional  problem  (projection  from  3D  to  2D),  then 
using  that  solution  to  set  up  a  one-dimensional  problem  in  the  transverse  direction  (projection  from 
3D  to  ID),  whose  solution  is  finally  used  to  reformulate  the  2D  minimization  problem. 

3.3.2  The  residuals  method 

The  residuals  method  consists  of  the  following  steps: 

Step  1:  Start  with  the  5-field  model  for  all  the  elements  in  the  mesh.  Obtain  the  finite  element  solution 
corresponding  to  this  model  as  described  in  Section  2.0.  This  step  is  the  same  as  in  the  previous 
method. 

Step  2:  Consider  a  candidate  model  2  from  the  available  set  in  Eq.  (40).  For  example,  we  may  con¬ 
sider  the  following  displacement  field: 

ur(r,s,w)  =  m|q(t,  ^)<J>1(w)  +  M|j(r,  5)03(w)  +  M|2(r,5)<l>6(w) 

us(r,s,w)  =  5)02(w)  +  «|j(r,  5)0>4(w)  +  m|2(m)07(w)  (42) 

uw(r,s,w)  =  u\™(r,  s)0>5(w)  +  m|^ (r,^)d>n(w) 

from  which  is  possible  to  compute  the  stiffness  matrix  and  load  vector  using  the  displacement  field  of 
Eq.  (42)  into  Eq.  (10).  The  resulting  system  of  equations  can  be  written  as: 

[^](2a){aJ"(2a)  =  ^{2  a)  (43) 

where  the  subscript  2 a  indicates  that  Eq.  (43)  refers  to  the  system  of  equations  for  candidate  model  2 
combination  a.  Note  that  the  size  of  the  solution  vector  {a}  is  larger  than  that  of  Eq.  (14)  because  the 
number  of  fields  is  now  eight  instead  of  five.  Compute  the  residuals  which  are  defined  as: 

'(2a)  =  {Q}(2a)~[K](2a){<*}l  (44) 

where  {a}\  is  the  solution  vector  corresponding  to  model  1  augmented  by  zeros  in  those  positions 
where  the  terms  corresponding  to  model  2  are  needed. 

Step  3:  Repeat  the  procedure  for  other  candidate  model  2,  in  such  a  way  that  the  residuals  for  each 
one  can  be  computed:  { R}{lb),  {/?}(2c),  iR}(2dy- Compute  the  norm  (vector  length)  of  each 
residual.  The  best  candidate  is  the  one  with  the  smallest  residual  norm. 
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Note  that  the  solution  of  the  system  of  equations  in  (43)  is  not  needed,  only  the  stiffness  matrices  and 
load  vectors  of  the  candidate  models  are  required. 

3.3.3  The  potential  energy  method 

This  criterion  for  model  selection  is  based  on  the  change  in  the  value  of  the  total  potential  energy  of 
the  problem.  The  combination  of  additional  fields  that  results  in  the  smallest  potential  energy  pro¬ 
vides  the  best  improvement  over  the  solution  of  model  1 . 

The  procedure  for  selecting  the  optimal  distribution  of  fields  for  a  given  problem  consists  of  two 
steps: 

Step  1:  The  basic  model  (5-field  semi-discretization)  is  increased  by  adding  one  more  field  to  each 
displacement  component  ( un  us,  uw),  in  such  a  way  that  the  transverse  function  for  the  r  and  s  compo¬ 
nents  have  the  same  order  of  director  function.  In  other  words,  model  2  should  be  an  8-field  semi-dis¬ 
cretization. 

Step  2:  The  criterion  for  selecting  which  of  the  available  director  functions  should  be  added  is  based 
on  the  change  in  the  value  of  the  total  potential  energy  of  the  problem.  Those  directors  functions 
which  result  in  the  smallest  potential  energy  will  provide  the  best  improvement  over  the  solution  of 
the  5-field  model.  The  potential  energy  accounts  for  the  effects  of  topology,  material  properties  and 
boundary  conditions,  thus  characterizing  the  problem.  The  potential  energy  is  computed  by  solving 
the  problem  several  times  for  each  candidate  combination  of  director  functions  at  a  low  p-level: 

n  =  \{af  [K]{a}-{af  [ q ] 


To  illustrate  the  concept  of  model  selection,  consider  the  situation  of  selecting  an  8-field  model.  There 
are  four  possible  combinations  from  the  expansion  given  in  Eq.  (40)  for  the  additional  three  fields 
needed  to  extend  the  5-field  model:  (a)  <D6,  <D7,  <b8;  (b)  06,  <E>7,  (c)  <D9,  O10,  d>8  and  (d)  <I>9,  <E>10, 

On.  No  other  combinations  are  possible,  given  the  constraint  on  the  order  of  the  director  functions 
indicated  before.  For  example,  the  hierarchic  model  2,  combination  (d)  would  be: 

ur(r,s,w )  =  u\rQ(r,  j)Oj(w)  +  u|J(r,  s)<f>3(w)  +  u\r2(r,s)<S*9(w) 
us(r,s,w )  =  w|J(r,  j)®2(w)  +  a|j(r,j)<&4(w)  +  i/|2(r,j)®10(w) 
uw(r,s,w)  =  u\™(r,  s)05(w)  +  wf^O^w) 

Of  all  the  methods  evaluated,  the  one  based  on  the  minimization  of  the  total  potential  energy  of  the 
problem  was  selected  for  implementation  into  the  prototype  software  because  it  was  the  simplest  and 
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most  effective  of  the  ones  analyzed.  The  computational  effort  to  determine  the  optimal  set  of  director 
functions  is  minimized  when  using  this  method.  The  automatic  procedure  implemented  in  the  proto¬ 
type  software  will  solve  for  each  combination  at  p-level=4,  and  select  that  one  which  minimizes  the 
potential  energy  of  the  problem.  Once  the  optimal  model  is  identified,  a  p-extension  is  performed  to 
ascertain  the  discretization  error. 

3.4  Example  problems 

3.4.1  Problem  4.  Thick  4-ply  laminate 

Consider  a  4-ply  [0/90]  s  laminated  plate  loaded  by  a  sinusoidal  transverse  load  qz(x,y)  =  -cos(nx/a) 
cos(ny/b).  The  four  layers  of  the  laminate  are  of  the  same  material  and  thickness  with  the  following 
properties: 

El  =  138000  MPa;  ET=  9300  MPa;  GLT  =  4600  MPa;  GTT=  3100  MPa;  vLT=  0.3,  v7T  =  0.5 

where  L  indicates  the  direction  parallel  to  the  fibers  and  T  is  the  transverse  direction.  When  the  L- 
direction  coincides  with  the  x-direction,  we  refer  to  it  as  the  0°  orientation  (Figure  12).  All  the  dimen¬ 
sions  are  in  millimeters. The  plate  is  hard-simply  supported  along  all  four  edges.  A  hard-simple  sup¬ 
port  is  characterized  by:  uz  =  ut  =  0,  where  ut  is  the  displacement  component  tangent  to  the  edge,  and 


^ - ► 


FIGURE  12.  Problem  4  -  Laminated  plate.  Notation 


The  reference  solution  was  obtained  from  the  finite  element  analysis  of  a  3D-solid  model  in  which 
each  layer  was  discretized  as  a  solid  element.  Because  of  symmetry,  only  one  quarter  of  the  plate  was 
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used  for  the  shell  analysis,  and  one  eighth  for  the  solid  analysis.  Figure  13  shows  the  one-element 


FIGURE  13.  Problem  4  -  Finite  element  meshes  for  the  shell  and  3D-solid 
models. 


mesh  for  the  shell  analysis  and  the  two-element  mesh  for  the  3D-solid  analysis.  Antisymmetry  con¬ 
straints  were  specified  on  the  middle  surface  of  the  3D-solid  model. 

The  results  for  the  shell  models  were  obtained  for  the  first  and  second  hierarchic  shell  models.  The 
optimal  combination  of  transverse  shape  functions  for  model  2  to  be  used  in  the  8-field  semi-discreti¬ 
zation  of  Eq.  (40)  was  determined  to  be  <J>6, 07,  On  (as  it  should  be  expected  for  a  bending  domi¬ 
nated  problem). 

The  results  shown  in  Table  3  include  the  total  potential  energy  of  the  solution  II(n);  the  maximum  dis¬ 
placement  at  the  center  of  the  plate  wz(0,0,0);  the  normal  stresses  ox(0,0,-h/2)  and  ay(0,0,-h/2)  and  the 
shear  stress  Txy(a/2,b/2,-h/2)  at  one  of  the  external  surfaces  of  the  plate  for  both  hierarchic  shell  models 
(shell  1  and  2)  and  for  the  3D-solid  model. 

table  3.  Results  for  problem  4 


Model 

n(tt)  x  io4 

Mz(0,0,0) 

Ox(0,0,-h/2) 

Cty(0,0,-h/2) 

Txy(a/2,b/2,-h/2) 

Shell  1 

-2.8746631 

-5.2565  x  10‘4 

6.4389 

1.4035 

-0.7274 

Shell  2 

-2.9863590 

-5.5671  x  10'4 

7.3924 

1.9688 

-0.8154 

3D-Solid 

-3.0850379 

-5.7263  x  10'4 

7.7388 

1.9267 

-0.8602 
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The  low  aspect  ratio  of  the  plate  combined  with  the  highly  anisotropic  nature  of  the  material  represent 
a  severe  test  for  any  laminated  shell  model.  The  length-to-thickness  ratio  is  only  3.5,  and  therefore  not 
suitable  for  conventional  shell  analysis.  However,  the  use  of  higher  order  models  indicates  that  the 
results  converge  to  the  fully  three-dimensional  solution. 

The  through-thickness  stress  distribution  of  the  normal  stresses  at  the  center  of  the  plate  is  shown  in 
Figure  14  for  shell  models  1  and  2.  Differences  between  the  two  models  is  visible  at  the  external  sur- 


FIGURE 14.  Problem  4  -  Through-thickness  stress  distributions  gx(0,0,.z)  and  oy(0,0,.z)  for  shell 
models  1  and  2. 

faces  and  at  the  interface  between  layers. 

3.4.2  Problem  5.  Effect  of  boundary  layers 

Boundary  layer  effects  occur  at  the  shell  boundaries,  and  are  characterized  by  the  fact  that  the  solution 
‘near’  the  boundary  is  substantially  different  from  the  solution  in  the  interior.  All  hierarchic  shell 
models  (as  well  as  the  fully  three-dimensional  model)  exhibit  boundary  layers,  and  an  important  part 
of  the  energy  of  the  solution  is  contained  in  them.  For  further  information  on  boundary  layer  effects 
refer  to  Refs.  [20]  to  [23].  Therefore,  the  mesh  design  necessary  to  obtain  accurate  solutions  for  any 
given  member  of  the  hierarchic  sequence  of  models  should  properly  account  for  the  boundary  layers. 
Extensive  numerical  experimentation  clearly  showed  that  the  hierarchic  models  are  very  capable  of 
resolving  the  boundary  layer  effects  when  proper  meshing  is  used. 


Design  and  Analysis  of  Composite  Multilayered  Shells 


35  of  52 


Higher  order  models 


Based  on  the  numerical  evidence,  guidelines  for  mesh  design  to  be  used  with  the  hierarchic  models 
that  will  provide  optimal  or  near-optimal  meshes  with  respect  to  the  energy  norm  were  developed. 
These  guidelines  are  summarized  as  follows: 

•  The  first  step  is  to  design  a  finite  element  mesh  that  provides  optimal  rate  of  convergence  for  the 
exact  solution  of  the  shell  problem  in  the  interior  of  the  domain  without  consideration  of  the 
boundaries.  For  smooth  problems  and  p-convergence  this  typically  involves  the  use  of  uniform  or 
quasi-uniform  meshes.  This  will  be  referred  to  as  the  “coarse  mesh”. 

•  Once  the  coarse  mesh  is  available,  the  boundary  layers  should  be  accounted  for  by  the  use  of 
graded  meshes.  For  most  practical  problems  one  or  two  layers  of  graded  elements  towards  the 
edges  are  sufficient  to  account  for  boundary  layer  effects.  The  characteristic  length  of  a  shell  prob¬ 
lem  is  the  thickness-to-radius  ratio  (h/R).  For  thin  shells  (h/R«  1),  the  recommended  size  of  the 
boundary  layer  elements  are  5  h  and  2>Jh  for  the  first  and  second  layers,  respectively.  For  thick 
shells  ( h/R  =  1),  boundary  layer  effects  are  less  significant  and,  in  general,  one  layer  of  elements  is 
sufficient  with  a  size  of  order  h. 

To  illustrate  the  effect  of  the  boundary  layer  in  the  solution  of  shell  problems,  consider  a  cylinder  with 
no  kinematical  constraints  at  the  ends,  subjected  to  a  sinusoidal  distributed  surface  traction 
Tn  =  rocos(20) .  This  load  is  self  equilibrated  in  the  angular  direction  and  uniform  in  the  axial 
direction.  The  dimensions  are  shown  in  Figure  15. 


FIGURE  15.  Problem  5  -  Cylinder  under  sinusoidal  loading.  Notation. 


Three  thickness-to-radius  ratios  were  analyzed:  h/R- 0.1,  h/R= 0.01,  and  h/R  =0.001  and  two  different 
materials  were  considered.  An  isotropic  material  with  £=10xl06,  v=l/3  and  unit  shear  factor;  and  a  4- 
ply  laminated  composite  with  a  [0/90]s  layup  and  the  following  properties  for  each  layer: 

El  =  25x106,  Et=  lxlO6,  Glt=  5x105,  G7T  =  2x105,  vlt=  0.25,  \TT=  0.49 
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where  L  indicates  the  direction  parallel  to  the  fibers  and  T  is  the  transverse  direction.  When  the  In¬ 
direction  coincides  with  the  Z-direction,  we  refer  to  it  as  the  0°  orientation:  For  the  two  outer  layers 
the  fibers  run  parallel  to  the  Z-axis,  and  for  the  two  inner  layers  the  fibers  run  in  the  circumferential 
direction.  All  layers  are  of  the  same  thickness  (7i/4).  This  problem  is  discussed  in  Ref.  [24]  for  the 
case  of  isotropic  material. 

The  radius  (R=1.0)  and  length  (L=2.0)  of  the  cylinder  are  kept  fixed,  and  the  thickness  is  changed  for 
each  case  analyzed.  The  solutions  were  obtained  for  hierarchic  models  1  and  2,  polynomial  orders 
ranging  from  1  to  8,  downward  run,  and  the  product  space  was  used.  For  definition  of  the  product 
space,  see  Ref.  [1],  page  96. 

Because  of  symmetry,  only  one  sixteenth  of  the  cylinder  is  considered  for  the  analysis.  Figure  16 
shows  the  finite  element  mesh  used  for  the  analysis  with  the  two  boundary  layer  elements  near  the 
free  end  of  the  shell  with  sizes  corresponding  to  h/R= 0.01,  that  is  £>]=(). 05  and  £>2=0-30. 


The  estimated  relative  error  in  energy  norm  as  a  function  of  the  number  of  degrees  of  freedom  (DOF) 
for  each  h/R  is  shown  in  Figure  17  for  the  case  of  isotropic  material  and  in  Figure  18  for  the  case  of 
the  4-ply  laminated  composite  material.  All  cases  shown  correspond  to  shell  model  1,  and  the  solu¬ 
tions  converged  to  less  than  1%  relative  error  in  energy  norm.  Note  that  the  rate  of  convergence  is 
very  low  (and  consequently  the  error  in  energy  norm  large)  for  polynomial  order  less  than  4  (run  #  5) 
for  the  thick  shell  (M?=0.1)  and  for  p-levels  less  than  5  (run  #4)  for  the  thin  shells  (h/R= 0.01  and 
0.001).  When  the  displacement  formulation  of  the  finite  element  method  is  used  for  thin  shells,  lock¬ 
ing  occurs  when  the  p-level  is  less  than  4  or  5,  depending  on  the  h/R  ratio. 
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Estimate,  33)=  SOL 

,  run  #1  to  #8 
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8 
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3 
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Convergence 
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2.31 

0.10 

1 

984 
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h/R  =  0.1 


■h/R=  0.01 


h/R  =  0.001 


FIGURE  17.  Problem  5  -  Estimated  relative  error  in  energy  norm.  Shell  model  1,  isotropic  case. 


The  boundary  layer  effects  can  be  visualized  when  displaying  the  first  principal  stress  distribution 
over  the  middle  surface  of  the  shell  (w  =  0.0).  Even  though  this  is  a  bending  dominated  problem,  the 
free-edge  boundary  layer  is  present  at  the  middle  surface  of  the  shell.  Figure  19  shows  the  first  princi¬ 
pal  stress,  S  j,  for  the  case  of  isotropic  material  and  for  all  three  thickness-to-radius  ratios.  Note  that  as 
the  thickness  of  the  shell  decreases,  Sj  is  practically  zero  everywhere,  except  along  a  narrow  band 
(the  boundary  layer)  near  the  free  edge,  and  that  the  size  of  the  boundary  layer  decreases  as  h/R  goes 
to  zero. 
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Error  Estimate,  U)=  SOL,  run  #1  to  #8 

Rate  of  Estimated 

Run  #  DOF  Total  Potential  Energy  Convergence  \  Error 

8  18  -5.314309984125694e-07  0.00  99.61 

7  66  -1.899446773714216e-05  0.12  84.97 

6  144  -6 . 512955636702421e-05  1.76  21.60 

S  252  -6 . 828094561142016e-05  4.02  2.28 

4  390  ~6 . 831639830084873e-05  6.31  0.15 

3  558  -6 . 831654203341285e-05  9.70  0.00 

2  756  -6 . 831654212361424e-05  1.75  0.00 

1  984  -6 . 831654215233987e-05  1.75  0.00 


Estimated  Limit 


-6 . 831654217133485e-05 


Error  Estimate,  XD=  SOL,  run  #1  to  #8 

Rate  of 

Run  #  DOF  Total  Potential  Energy  Convergence 


Estimated 
%  Error 


h/R  =0.1 


8 

18 

-5 . 367603163346339e-06 

0.00 

100.00 

7 

66 

-2 . 683182682004668e-04 

0.00 

99.79 

6 

144 

-4 . 745771588527027e-02 

0.83 

52.29 

5 

252 

-6 . 467322563428622e-02 

2.97 

9.94 

4 

390 

-6 . 531366484827916e-02 

5.56 

0.88 

3 

558 

-6 . 531869074802141e-02 

8.83 

0.04 

2 

756 

-6 . 531869475402023e-02 

0.97 

0.03 

1 

984 

-6 . 531869675536421e-02 

0.97 

0.02 

Estimated  Limit 

-6 . 531869976328750e-02 

Error 

Estimate,  ED= 

SOL,  run  #1  to  #8 

Rate  of 

Estimated 

Run  # 

DOF 

Total  Potential  Energy 

Convergence 

%  Error 

8 

18 

-5 . 368169592016818e-05 

0.00 

100.00 

7 

66 

-2 . 695259874825106e-03 

0.00 

100.00 

6 

144 

-4 . 590992118544692e+01 

0.78 

54.48 

5 

252 

-6 . 429500566500084e+01 

2.66 

12.32 

4 

390 

-6 . 526508159525335e+01 

4.44 

1.77 

3 

558 

-6 . 528518129525122e+01 

5.37 

0.26 

2 

756 

-6 . 528560184167117e+01 

5.28 

0.05 

1 

984 

-6 . 528561843952617e+01 

5.28 

0.01 

Estimated  Limit 

-6 . 528561953032957e+01 

FIGURE  18.  Problem  5  -  Estimated  relative  error  in  energy  norm.  Shell  model  1 ,  laminated  composite 
case. 
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FIGURE  19.  Problem  5  -  First  principal  stress  distribution  at  the  middle  surface  of  the  shell.  Shell  model  1 , 
isotropic  case. 
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The  boundary  layer  effect  is  also  present  on  the  external  surface  ( w  =  hi 2)  of  the  shell.  Considering 
the  case  of  isotropic  material  shown  in  Figure  20,  the  Sj  stress  distribution  is  rather  regular  every- 


i»*m-  i 

Mm*M-**?**W 

■ 

!**«•«»** 

IJiftfetM 

- 

a 

a 

: 

:a 

is&avfe 

'mm*  mvmm 
\wh* 


I 


LLJ 


ftjMfattft 


I  iMsunm 


FIGURE  20.  Problem  5  -  First  principal  stress  distribution  at  w=h/2.  Shell  model  1,  isotropic  case. 


where  except  close  to  the  shell  free  end.  where  the  presence  of  the  boundary  layer  perturb  the  stress 
distribution. 

The  situation  is  quite  similar  for  the  case  of  the  laminated  composite  shell.  Figure  21  shows  the  Sj 
stress  distribution  for  the  middle  surface  of  the  shell,  that  is  at  w  =  0.0,  for  all  three  h/R  ratios.  Com¬ 
paring  Figure  19  with  Figure  21,  the  behavior  of  the  boundary  layer  is  almost  identical  for  both  mate¬ 
rials. 
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FIGURE  21.  Problem  5  -  First  principal  stress  distribution  at  the  middle  surface  of  the  shell.  Shell  model  1 , 
Laminated  composite  case. 


The  effect  of  the  boundary  layer  is  less  apparent  at  the  interface  between  layers,  however.  As  shown 
in  Figure  22,  the  stress  distribution  near  the  free  edge  is  only  mildly  perturbed  by  the  boundary  layer. 
The  laminated  composite  case  was  also  analyzed  using  shell  model  2.  No  substantial  difference  can 
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FIGURE  22.  Problem  5  -  First  principal  stress  distribution  at  w-h/4.  Shell  model  1 ,  Laminated  composite 
case. 


be  realized  between  the  results  of  model  1  and  model  2,  however.  The  through-thickness  normal 
stress  distribution  Sy  for  h/R- 0.1,  at  a  point  located  at  0=0  on  the  free  end  of  the  cylinder  (point  A  in 
Figure  15),  is  shown  in  Figure  23.  The  results  for  shell  models  1  and  2  are  almost  indistinguishable 
from  each  other.  The  same  situation  was  found  for  the  other  two  h/R  ratios. 

The  first  principal  stress  distribution,  Sj,  at  the  middle  surface  of  the  laminated  composite  shell  is 
shown  in  Figure  24  for  h/R= 0.1  and  h/R= 0.01  as  computed  from  shell  model  2.  Comparing  this  plot 
with  those  shown  in  Figure  21,  it  is  clear  that  the  same  type  of  boundary  layer  is  present  in  both  shell 
models,  and  that  the  localized  stresses  induced  are  of  the  same  order  of  magnitude. 
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FIGURE  23.  Problem  5  -  Through-thickness  normal  stress  distribution  at  point  A.  Shell  models  1  and  2, 

laminated  composite  case. 

Table  4  shows  the  normalized  displacement  of  point  A  for  all  three  h/R  ratios  and  for  the  isotropic  and 
laminated  composite  cases.  The  4-ply  laminate  results  are  shown  for  shell  models  1  and  2.  The  nor¬ 
malized  displacement  is  defined  as: 

Ertfux 

U  =  — — / 

T0R 

where  Ej=  lxlO6,  T0=1.0,  R=l.O  and  uxA  is  the  displacement  of  point  A  in  the  global  x-direction 
(Figure  15).  Note  that  U  converges  to  a  limit  value  as  the  thickness-to  radius  ratio  goes  to  zero,  and 

table  4.  Problem  5  -  Normalized  displacement  of  point  A. 


Normalized  displacement  U 


h/R 

Isotropic  model  1 

Laminated  model  1 

Laminated  model  2 

0.1 

0.127 

0.349 

0.346 

0.01 

0.120 

0.333 

0.331 

0.001 

0.119 

0.333 

0.330 

that  both  hierarchic  models  converge  to  the  same  value. 


Design  and  Analysis  of  Composite  Multilayered  Shells 


44  of  52 


Higher  order  models 


&*VHb 

ft 

Mm  * 

life*  lilPISs*# 

■ 

!  MMkm' 

iff! 

in 

>!. 

■ 

l 

l 

■ 

■ 

1 

!  fiJOjMMHI 

mm*  wmmm 


jtJMtttKS 


wimm 


wvmm 

tit): 

IWI 


FIGURE  24.  Problem  5  -  First  principal  stress  distribution  at  the  middle  surface  of  the  shell.  Shell  model  2, 

Laminated  composite  case. 


3.4.3  Problem  6.  Effect  of  mesh  distortion 

The  same  model  problem  5  is  used  to  demonstrate  the  robustness  of  the  shell  models  in  the  presence 
of  distorted  meshes.  In  particular,  we  are  interested  in  evaluating  the  influence  of  the  element  edges 
not  being  aligned  with  the  principal  directions  of  the  shell  surface. 

Consider  the  cylindrical  shell  shown  in  Figure  15  with  the  finite  element  mesh  of  Figure  25.  The  6- 
element  mesh  was  designed  in  such  a  way  that  the  center  longitudinal  line  can  be  rotated  through  an 
arbitrary  angle  a  to  change  the  distortion  of  the  elements.  The  analysis  was  performed  for  the  isotro¬ 
pic  case  and  for  two  thickness-to-radius  ratios:  h/R= 0.01  and  h/R= 0.001.  Downward  p  extension  was 
used  together  with  the  product  space  and  shell  model  1 . 
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FIGURE  25.  Problem  6  -  Distorted  finite  element  mesh,  6  elements 


The  results  of  the  analysis  are  shown  in  Figure  26,  where  the  estimated  relative  error  in  energy  norm 
is  shown  for  three  values  of  the  distortion  angle  a.  Note  that  the  potential  energy  of  the  solution  is 
practically  independent  on  the  distortion  angle,  and  for  all  three  values  of  a  the  relative  error  in 
energy  norm  at  p-level=8  (run  #1)  is  very  small.  The  results  are  summarized  in  Table  5  for  both  h/R 
ratios.  Included  in  the  table  are  the  values  of  the  potential  energy  corresponding  to  a  p-level  of  8  (run 
#1)  and  the  ux  displacement  component  of  point  A  (see  Figure  15). 

table  5.  Problem  6  -  Effect  of  distortion  angle,  6-element  mesh 


h/R 

a  [deg] 

Potential  Energy 

0 

-2.399350588x1  O'2 

0.122194 

0.01 

5 

-2.39935 1950x1  O'2 

0.122194 

10 

-2.3993543 19xl0’2 

0.122194 

0 

-23.87654257 

121.602 

0.001 

5 

-23.87655967 

121.602 

10 

-23.87657944 

121.602 
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Error 

Estimate,  ID=  SOL,  run  #1  to  #8 

Rate  of 

Estimated 

Run  # 

DOF  Total  Potential  Energy 

Convergence 

\  Error 

8 

36  -1 . 336241124934671e-05 

0.00 

99.97 

7 

132  -2 . 84061S710935025e-03 

0.05 

93.89 

6 

288  -2 . 343051045818392e-02 

2.32 

15.32 

5 

504  -2 . 399169471407069e-02 

5.13 

0.87 

4 

780  -2 . 399347533110004e-02 

4.67 

0.11 

3 

1116  -2 . 399349983960644e-02 

2.25 

0.05 

2 

1512  -2 . 399350534096622e-02 

3.75 

0.02 

1 

1968  -2 . 399350588291860e-02 

3.75 

0.01 

Estimated  Limit  -2.399350597013211e-02 

Error 

Estimate,  ID^  SOL,  run  #1  to  #8 

Rate  of 

Estimated 

Run  # 

DOF  Total  Potential  Energy 

Convergence 

%  Error 

8 

36  -1 . 289755202986731e-05 

0.00 

99.97 

7 

132  -2 . 654009335612438e-03 

0.04 

94.31 

6 

288  -2 . 328926272443428e-02 

2.19 

17.13 

5 

504  -2 . 399084600980942e-02 

4.98 

1.06 

4 

780  -2 . 399348754963635e-02 

5.06 

0.12 

3 

1116  -2 . 399351370342892e-02 

2.36 

0.05 

2 

1512  -2 . 399351895602165e-02 

3.65 

0.02 

1 

1968  -2 . 399351950582236e-02 

3.65 

0.01 

Estimated  Limit  -2 . 399351960005294e-02 

Error 

Estimate,  ID-  SOL,  run  #1  to  #8 

Rate  of 

Estimated 

Run  # 

DOF  Total  Potential  Energy 

Convergence 

%  Error 

8 

36  -1 . 161274162268969e-05 

0.00 

99.98 

7 

132  -2 . 139050734588481e-03 

0.04 

95.44 

6 

288  -2 . 285085355920161e-02 

1.89 

21.82 

5 

504  -2 . 398729137464645e-02 

4.65 

1.61 

4 

780  -2 . 399350361420738e-02 

5.79 

0.13 

3 

1116  -2 . 399353761288973e-02 

2.71 

0.05 

2 

1512  -2 . 399354267072279e-02 

3.65 

0.02 

1 

1968  -2 . 399354319955519e-02 

3.65 

0.01 

Estimated  Limit  -2 . 399354329008791e-02 

oc  =  0° 


a=  5° 


a=  10° 


FIGURE  26.  Problem  6  -  Estimated  relative  error  in  energy  norm  for  h/R  =0.01 . 
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A  second  mesh  design  was  also  considered  for  this  problem  as  shown  in  Figure  27,  where  a  total  of  10 


FIGURE  27.  Problem  6  -  Distorted  finite  element  mesh,  10  elements. 


shell  elements  were  used.  Both  free  ends  of  the  cylinder  are  included  in  this  mesh,  and  therefore 
boundary  layer  elements  are  required  at  each  end. 

The  results  for  this  mesh  are  summarized  in  Table  6  for  h/R= 0.005. 

table  6.  Problem  6  -  Effect  of  distortion  angle,  10-element  mesh 


h/R 

a  [deg] 

Potential  Energy 

uxA 

0 

-0.38308005823 

0.976851 

0.005 

5 

-0.38308003711 

0.976851 

10 

-0.38308000626 

0.976852 

All  the  results  indicate  the  very  low  sensitivity  of  the  solution  to  element  distortion,  even  in  the  pres¬ 
ence  of  boundary  layers. 


Design  and  Analysis  of  Composite  Multilayered  Shells 


48  of  52 


Summary  and  Conclusions 


4.0  Summary  and  Conclusions 


All  objectives  indicated  in  the  Phase  II  proposal  have  been  achieved.The  results  of  the  research  effort 

were  implemented  in  the  prototype  software  which  utilizes  the  data  structure  of  Stress  Check.  The 

main  accomplishments,  described  in  detail  in  this  report,  can  be  summarized  as  follows: 

•  A  hierarchic  sequence  of  models  for  laminated  composite  shells  was  developed  and  implemented 
in  the  prototype  software  for  the  solution  of  linear  elastostatics  problems.  These  models  are  also 
capable  of  solving  problems  with  isotropic  or  orthotropic  homogeneous  materials.  Only  three  hier¬ 
archic  models  were  implemented  for  laminated  composites  during  this  research  Phase  II.  The  pro¬ 
totype  software  was  set  up  to  accommodate  more  models  once  they  become  available  in  the  final 
commercial  implementation  during  the  Phase  III  project. 

•  An  automatic  procedure  for  the  selection  of  higher  order  models  from  the  hierarchic  family  of 
models  was  developed  and  implemented  in  the  prototype  software.  The  automatic  procedure, 
based  on  the  change  in  the  total  potential  energy  of  the  problem,  can  be  disabled  so  that  the  models 
can  be  manually  selected  instead. 

•  A  unique  mapping  technique  for  shells  (quasi-regional  mapping)  was  implemented  in  the  proto¬ 
type  software  to  handle  completely  general  surface  descriptions.  The  surfaces  that  can  be  created 
in  the  prototype  software  include  spheres,  cylinders,  cones,  torus,  general  surfaces  created  by  using 
3D-curves:  Tabulated  cylinders,  extruded  surfaces,  ruled  surfaces  and  surface  of  revolution;  and 
surfaces  created  by  using  control  points:  Non-uniform  rational  B-  and  P-splines. 

•  The  research  work  leading  to  the  geometric  nonlinear  analysis  of  shells  was  completed  in  the  pro¬ 
totype  software,  and  the  concept  was  implemented  and  tested  in  the  3D-solids  environment.  The 
implementation  for  shells  will  be  addressed  during  the  Phase  HI  project.  Additionally,  the  research 
and  implementation  of  eigenvalue  buckling  and  prestress  modal  analyses  were  completed  during 
the  Phase  II  project. 

•  The  procedures  for  the  computation  of  stiffness  matrices,  mass  matrices,  geometric  matrices  and 
load  vectors  were  developed  and  implemented  in  the  prototype  software  for  quadrilateral  shell  ele¬ 
ments.  Similar  procedures  for  triangular  shell  elements  were  initiated  but  not  completed  during  the 
Phase  II  project.  These  are  implementation  issues,  not  research  issues,  and  therefore  postponed  for 
the  Phase  HI  project. 

•  The  following  boundary  conditions  were  implemented  and  tested  in  the  prototype  software  for  all 
the  members  of  the  hierarchic  sequence  of  models: 

Loads:  Surface  tractions  can  be  specified  in  global  or  local  coordinate  system  or  in  the  direction 
normal  to  the  shell  surface.  Edge  tractions  can  be  specified  as  normal/tangent  membrane,  bend¬ 
ing/twisting  moments  and  transverse  shear.  Body  forces  can  also  be  assigned  in  the  global  coor¬ 
dinate  directions.  All  applied  loads  can  be  defined  as  constant,  parametric  or  formula.  Point 
loads  can  be  applied  to  a  node,  point,  or  in  the  interior  of  the  shell  element  either  in  the  global  or 
in  a  local  reference  frame.  The  polynomial  order  of  the  elements  with  points  loads  should  be 
kept  low  (typically  p=3  or  4). 
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Constraints:  Edge  constraints  can  be  specified  as  simply  support,  symmetry,  antisymmetry, 
sliding  support,  pinned  support  or  built-in.  Node  constraints  are  allowed  to  prevent  rigid  body 
displacement  and  rotation  only. 

•  Several  representative  models  problems  were  solved  using  the  hierarchic  models,  and  their  results 
compared  with  those  available  in  the  literature  or  with  reference  solutions  obtained  by  solving  3D- 
solid  models.  The  results  clearly  show  the  ability  of  the  shell  models  to  approximate  the  three- 
dimensional  problem  well.  Many  more  problems  than  the  ones  included  in  this  report  were  ana¬ 
lyzed,  but  only  a  representative  set  was  selected  to  illustrate  the  main  points  of  the  implementation. 

The  phase  II  project  clearly  demonstrated  the  potential  of  the  use  of  a  hierarchic  sequence  of  models 
for  the  analysis  and  design  of  composite  multilayered  shells.  Phase  III  project  will  address  the  com¬ 
mercial  implementation  of  these  models  so  that  the  engineering  community  will  benefit  from  this  new 
technology. 
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